# Deep Unfolded Tensor Robust PCA with Self-supervised Learning

Harry Dong
CMU
Department of Electrical and Computer Engineering, Carnegie Mellon University,
Pittsburgh, PA 15213, USA; Emails<EMAIL_ADDRESS> Megna
Shah
AFRL
Materials and Manufacturing Directorate, Air Force Research Laboratory,
Wright-Patterson AFB, OH 45433, USA; Emails:
<EMAIL_ADDRESS>Sean Donegan22footnotemark: 2
AFRL
Yuejie Chi11footnotemark: 1
CMU

###### Abstract

Tensor robust principal component analysis (RPCA), which seeks to separate a
low-rank tensor from its sparse corruptions, has been crucial in data science
and machine learning where tensor structures are becoming more prevalent.
While powerful, existing tensor RPCA algorithms can be difficult to use in
practice, as their performance can be sensitive to the choice of additional
hyperparameters, which are not straightforward to tune. In this paper, we
describe a fast and simple self-supervised model for tensor RPCA using deep
unfolding by only learning four hyperparameters. Despite its simplicity, our
model expunges the need for ground truth labels while maintaining competitive
or even greater performance compared to supervised deep unfolding.
Furthermore, our model is capable of operating in extreme data-starved
scenarios. We demonstrate these claims on a mix of synthetic data and real-
world tasks, comparing performance against previously studied supervised deep
unfolding methods and Bayesian optimization baselines.

Keywords: tensors, robust principal component analysis, self-supervised
learning, deep unfolding.

##  1 Introduction

Becoming increasingly pervasive in data science and machine learning, tensors
are powerful data structures that capture interactions between elements in
higher dimensions. These structures lead to inherent properties that can be
exploited for tasks that seek to analyze and process tensors. One such task is
robust principal component analysis (RPCA), the problem of recovering a low-
rank tensor that has been sparsely corrupted, which has found numerous
applications such as surveillance, anomaly detection, and more.

Compared with its matrix counterpart, dealing with tensor RPCA faces some
unique challenges, as many ideas from matrix RPCA begin to fall apart. First,
blindly applying a matrix RPCA algorithm to a flattened tensor can destroy
structural information [YZ16]. Second, convex formulations, typically done
with the nuclear norm as a convex relaxation to the rank constraint, is NP-
hard for tensors [FL18]. Third, there are multiple notions of tensor
decomposition and rank. As such, provable tensor RPCA algorithms have been
developed for low multilinear rank [DTMC22], tubal rank [LFC+16, LFC+19], and
CP-rank [DBBG19, AJSN16] decompositions.

In this paper, we focus on the Tucker decomposition with low multilinear rank,
though we believe our approach is also extendable to other tensor
decompositions. While powerful, existing efficient tensor RPCA algorithms,
e.g. based on iterative scaled gradient updates as in [DTMC22], can be
difficult to use in practice, as their performance is sensitive to the choice
of hyperparameters, such as learning rates and thresholds. In addition, tuning
these hyperparameters can be difficult, since the number of hyperparameters
scales with the desired number of iterations.

Leveraging deep unfolding, the process of unrolling an iterative algorithm
into a deep neural network introduced in [GL10], there have been considerable
efforts in learning algorithmic parameters with backpropagation to improve the
performance of the original algorithm. Deep unfolding for RPCA has been
applied to ultrasound imaging [SCZ+19, CZS+19], background subtraction
[VLJED21], and the special case of positive semidefinite (PSD) low-rank
matrices [HKK+20]. To generalize this idea, [CLY21] designed a deep unfolded
architecture for matrix RPCA, which is extendable to infinitely many RPCA
iterations. However, many existing approaches [SCZ+19, CZS+19, VLJED21, CLY21]
train on ground truth labels which may be limited or nonexistent in practice.
Though only applied to PSD low-rank matrices, [HKK+20] addressed this issue
with an unsupervised model.

Motivated by the need to broaden the applicability and improve the performance
of tensor RPCA in practice, we describe a fast and simple self-supervised
model for tensor RPCA by unfolding the recently proposed algorithm in
[DTMC22], due to its appealing performance. Our contributions are as follows.

  1. 1.

We propose a novel self-supervised model for tensor RPCA, which can be used
independently or as fine tuning for its supervised counterpart. Furthermore,
it scales to an arbitrary number of RPCA iterative updates with only four
learnable parameters by leveraging theoretical insights from [DTMC22].

  2. 2.

Synthetic and real-world experiments on video surveillance and materials
microscopy data show that our self-supervised model matches or exceeds the
performance of supervised learning without the need for ground truth,
especially in data-starved scenarios.

#### Paper organization.

The rest of this paper is organized as follows. We build up the tensor RPCA
algorithm from [DTMC22] in Section 2 which will be used in our learned tensor
RPCA method described in Section 3. Then, we present synthetic and real-world
experimental results in Section 4, followed by final remarks in Section 5.

#### Notation and basics of tensor algebra.

Throughout this paper, we represent tensors with bold calligraphic letters
(e.g. $\bm{\mathcal{X}}$) and matrices with bold capitalize letters (e.g.
$\bm{X}$). For tensor $\bm{\mathcal{A}}\in\mathbb{R}^{n_{1}\times n_{2}\times
n_{3}}$, $\mathcal{M}_{k}\left(\bm{\mathcal{A}}\right)$ is the tensor
matricization along the $k$-th mode (i.e. dimension of the tensor), $k=1,2,3$,
and $[\bm{\mathcal{A}}]_{i_{1},i_{2},i_{3}}$ is the $(i_{1},i_{2},i_{3})$-th
entry of $\bm{\mathcal{A}}$. A fiber is a vector obtained by fixing all
indices except one in the tensor (e.g.
$[\bm{\mathcal{A}}]_{i_{1},\cdot,i_{3}}$). Let the inner product between two
tensors $\bm{\mathcal{A}}$ and $\bm{\mathcal{B}}$ be
$\left\langle\bm{\mathcal{A}},\bm{\mathcal{B}}\right\rangle=\sum_{i_{1},i_{2},i_{3}}[\bm{\mathcal{A}}]_{i_{1},i_{2},i_{3}}[\bm{\mathcal{B}}]_{i_{1},i_{2},i_{3}}$,
and
$\left\lVert\bm{\mathcal{A}}\right\rVert_{\operatorname{\mathsf{F}}}=\sqrt{\left\langle\bm{\mathcal{A}},\bm{\mathcal{A}}\right\rangle}$
and
$\left\lVert\bm{\mathcal{A}}\right\rVert_{1}=\sum_{i_{1},i_{2},i_{3}}|[\bm{\mathcal{A}}]_{i_{1},i_{2},i_{3}}|$
denote the Frobenius norm and the $\ell_{1}$-norm, respectively. Furthermore,
suppose a tensor $\bm{\mathcal{X}}\in\mathbb{R}^{n_{1}\times n_{2}\times
n_{3}}$ has multilinear rank $\bm{r}=(r_{1},r_{2},r_{3})$, and its Tucker
decomposition is given by
$\bm{\mathcal{X}}=(\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)})\operatorname{\bm{\cdot}}\bm{\mathcal{G}}$
with $\bm{U}^{(k)}\in\mathbb{R}^{n_{k}\times r_{k}}$ and
$\bm{\mathcal{G}}\in\mathbb{R}^{r_{1}\times r_{2}\times r_{3}}$, it follows
that

$\displaystyle[\bm{\mathcal{X}}]_{i_{1},i_{2},i_{3}}=\sum_{j_{1},j_{2},j_{3}}\Bigg{(}\prod_{k=1,2,3}[\bm{U}^{(k)}]_{i_{k},j_{k}}\Bigg{)}[\bm{\mathcal{G}}]_{j_{1},j_{2},j_{3}}.$

Last but not least, given any $\bm{\mathcal{X}}\in\mathbb{R}^{n_{1}\times
n_{2}\times n_{3}}$, its rank-$\bm{r}$ higher-order singular value
decomposition (HOSVD) $\mathcal{H}_{\bm{r}}\left(\bm{\mathcal{X}}\right)$ is
given by
$\mathcal{H}_{\bm{r}}\left(\bm{\mathcal{X}}\right)=(\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)},\bm{\mathcal{G}})$,
where $\bm{U}^{(k)}\in\mathbb{R}^{n_{k}\times r_{k}}$ contains the top $r_{k}$
singular vectors of $\mathcal{M}_{k}(\bm{\mathcal{X}})$ and
$\bm{\mathcal{G}}=(\bm{U}^{(1)\top},\bm{U}^{(2)\top},\bm{U}^{(3)\top})\operatorname{\bm{\cdot}}\bm{\mathcal{X}}$.

##  2 Background on Tensor RPCA

###  2.1 Problem formulation

Suppose we observe a corrupted tensor
$\bm{\mathcal{Y}}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ of the form

$\displaystyle\bm{\mathcal{Y}}=\bm{\mathcal{X}}_{\star}+\bm{\mathcal{S}}_{\star},$
(1)

where $\bm{\mathcal{X}}_{\star}$ is a low-rank tensor with the multilinear
rank $\bm{r}=(r_{1},r_{2},r_{3})$, whose Tucker decomposition is given by
$\bm{\mathcal{X}}_{\star}=(\bm{U}_{\star}^{(1)},\bm{U}_{\star}^{(2)},\bm{U}_{\star}^{(3)})\operatorname{\bm{\cdot}}\bm{\mathcal{G}}_{\star}$
with $\bm{U}_{\star}^{(k)}\in\mathbb{R}^{n_{k}\times r_{k}}$ for $k=1,2,3$,
and $\bm{\mathcal{G}}_{\star}\in\mathbb{R}^{r_{1}\times r_{2}\times r_{3}}$.
Moreover, $\bm{\mathcal{S}}_{\star}$ is an $\alpha$-sparse corruption tensor,
i.e. $\bm{\mathcal{S}}_{\star}$ contains at most $0\leq\alpha<1$ fraction
nonzero values along each fiber. Given $\bm{\mathcal{Y}}$ and $\bm{r}$, the
goal of tensor RPCA aims to accurately obtain $\bm{\mathcal{X}}_{\star}$ and
$\bm{\mathcal{S}}_{\star}$.

###  2.2 Tensor RPCA via ScaledGD

Recently, Dong et al. [DTMC22] proposed a fast and scalable method for tensor
RPCA with provable performance guarantees. Specifically, they try to minimize
the objective function

$\displaystyle\mathcal{L}(\bm{F},\bm{\mathcal{S}})\coloneqq\frac{1}{2}\left\|\big{(}\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)}\big{)}\operatorname{\bm{\cdot}}\bm{\mathcal{G}}+\bm{\mathcal{S}}-\bm{\mathcal{Y}}\right\|_{\operatorname{\mathsf{F}}}^{2},$
(2)

where $\bm{F}=(\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)},\bm{\mathcal{G}})$ and
$\bm{\mathcal{S}}$ are the estimates of the factors of the low-rank tensor and
the sparse tensor, respectively. The algorithm proceeds by iterative updates
of the factors $\bm{F}$ via scaled gradient descent (ScaledGD) [TMC21a,
TMPB+22, TMC22, TMC21b], followed by filtering the outliers $\bm{\mathcal{S}}$
using the shrinkage operator $\mathcal{T}_{\zeta}(\cdot)$, defined as

$\displaystyle[\mathcal{T}_{\zeta}(\bm{\mathcal{X}})]_{i,j,k}\coloneqq\mathrm{sgn}([\bm{\mathcal{X}}]_{i,j,k})\cdot\max(0,|[\bm{\mathcal{X}}]_{i,j,k}|-\zeta).$
(3)

The details of the proposed ScaledGD algorithm are summarized in Algorithm 1.
Beginning with a careful initialization using the spectral method, given by

$\displaystyle\big{(}\bm{U}_{0}^{(1)},\bm{U}_{0}^{(2)},\bm{U}_{0}^{(3)},\bm{\mathcal{G}}_{0}\big{)}=\mathcal{H}_{\bm{r}}\left(\bm{\mathcal{Y}}-\mathcal{T}_{\zeta_{0}}\left(\bm{\mathcal{Y}}\right)\right),$
(4)

for each iteration, ScaledGD makes the following updates:

$\displaystyle\bm{\mathcal{S}}_{t+1}$
$\displaystyle=\mathcal{T}_{\zeta_{t+1}}\left(\bm{\mathcal{Y}}-\big{(}\bm{U}_{t}^{(1)},\bm{U}_{t}^{(2)},\bm{U}_{t}^{(3)}\big{)}\operatorname{\bm{\cdot}}\bm{\mathcal{G}}_{t}\right),$
(5a) $\displaystyle\bm{U}^{(k)}_{t+1}$
$\displaystyle=\bm{U}^{(k)}_{t}-\eta\nabla_{\bm{U}^{(k)}_{t}}\mathcal{L}(\bm{F}_{t},\bm{\mathcal{S}}_{t+1})\big{(}\breve{\bm{U}}_{t}^{(k)\top}\breve{\bm{U}}_{t}^{(k)}\big{)}^{-1}$
(5b) for $k=1,2,3$, and $\displaystyle\bm{\mathcal{G}}_{t+1}$
$\displaystyle=\bm{\mathcal{G}}_{t}-\eta\left(\tilde{\bm{U}}^{(1)}_{t},\tilde{\bm{U}}^{(2)}_{t},\tilde{\bm{U}}^{(3)}_{t}\right)\operatorname{\bm{\cdot}}\nabla_{\bm{\mathcal{G}}_{t}}\mathcal{L}(\bm{F}_{t},\bm{\mathcal{S}}_{t+1}),$
(5c)

where $\eta>0$ is the learning rate, and $\\{\zeta_{t}\\}_{t=0}^{T}$ is the
threshold schedule up to $T$ iterations. Here,

$\displaystyle\breve{\bm{U}}_{t}^{(1)}$
$\displaystyle=(\bm{U}_{t}^{(3)}\otimes\bm{U}_{t}^{(2)})\mathcal{M}_{1}(\bm{\mathcal{G}}_{t})^{\top},$
$\displaystyle\breve{\bm{U}}_{t}^{(2)}$
$\displaystyle=(\bm{U}_{t}^{(3)}\otimes\bm{U}_{t}^{(1)})\mathcal{M}_{2}(\bm{\mathcal{G}}_{t})^{\top},$
$\displaystyle\breve{\bm{U}}_{t}^{(3)}$
$\displaystyle=(\bm{U}_{t}^{(2)}\otimes\bm{U}_{t}^{(1)})\mathcal{M}_{3}(\bm{\mathcal{G}}_{t})^{\top},$

and
$\tilde{\bm{U}}^{(k)}_{t}=\big{(}\bm{U}^{(k)\top}_{t}\bm{U}_{t}^{(k)}\big{)}^{-1}$,
$k=1,2,3$, where $\otimes$ denotes the Kronecker product. In [DTMC22], it is
demonstrated that ScaledGD allows perfect recovery of the low-rank tensor
$\bm{\mathcal{X}}_{\star}$ and the sparse tensor $\bm{\mathcal{S}}_{\star}$
under mild assumptions as long as the low-rank tensor is incoherent and the
corruption level $\alpha$ is not too large. Moreover, ScaledGD converges at a
linear rate independent of the condition number of $\bm{\mathcal{X}}_{\star}$
[DTMC22], making it more appealing than the vanilla gradient descent approach
that is much more sensitive to ill-conditioning. However, translating the
theoretical advantage into practice requires carefully tuned hyperparamters,
$\eta$ and $\\{\zeta_{t}\\}_{t=0}^{T}$, which unfortunately are not readily
available.

Algorithm 1 ScaledGD for 3rd order tensor RPCA

Input: the observed tensor $\bm{\mathcal{Y}}$, the multilinear rank $\bm{r}$,
learning rate $\eta$, and threshold schedule $\\{\zeta_{t}\\}_{t=0}^{T}$.

Init:
$\big{(}\bm{U}_{0}^{(1)},\bm{U}_{0}^{(2)},\bm{U}_{0}^{(3)},\bm{\mathcal{G}}_{0}\big{)}=\mathcal{H}_{\bm{r}}\left(\bm{\mathcal{Y}}-\mathcal{T}_{\zeta_{0}}\left(\bm{\mathcal{Y}}\right)\right)$.

for $t=0,1,\dots,T-1$ do

Update $\bm{\mathcal{S}}_{t+1}$ via (5a);

Update
$\bm{F}_{t+1}=\big{(}\bm{U}^{(1)}_{t+1},\bm{U}^{(2)}_{t+1},\bm{U}^{(3)}_{t+1},\bm{\mathcal{G}}_{t+1}\big{)}$
via (5b) and (5c);

end for

Output:
$\bm{F}_{T}=\big{(}\bm{U}_{T}^{(1)},\bm{U}_{T}^{(2)},\bm{U}_{T}^{(3)},\bm{\mathcal{G}}_{T}\big{)}$.

##  3 Proposed Method

We aim to greatly broaden the applicability of Algorithm 1 with learned
hyperparameters by leveraging self-supervised learning (SSL) under a deep-
unfolding perspective [GL10] of Algorithm 1.

###  3.1 Unrolling-aided hyperparameter tuning

Based on the insights from theoretical analysis [DTMC22, Theorem 1], we use a
threshold schedule that decays at every iteration to reduce the number of
hyperparameters. In other words,

$\zeta_{t+1}=\rho\zeta_{t}$

for $t\geq 1$ and $0<\rho<1$, so there are only 4 hyperparameters that need
tuning: $\zeta_{0}$, $\zeta_{1}$, $\rho$, and $\eta$. By unrolling Algorithm
1, we obtain an unfolded version of Algorithm 1, where $\zeta_{0}$,
$\zeta_{1}$, $\rho$, and $\eta$ are learnable; see Fig. 1 for more details.
The model consists of a feed-forward layer and a recurrent layer to model the
initialization and iterative steps of Algorithm 1.

Figure 1: Network architecture. $\bm{\mathcal{Y}}$ and $\bm{r}$ are fed into
the first layer, the initialization step in Algorithm 1. Subsequent recurrent
layers follow Algorithm 1 to update
$\bm{\mathcal{X}}_{t}=(\bm{U}_{t}^{(1)},\bm{U}_{t}^{(2)},\bm{U}_{t}^{(3)})\operatorname{\bm{\cdot}}\bm{\mathcal{G}}_{t}$.
Recurrent layers require the iteration count $t$ due to our definition of the
$(t+1)$-th thresholding parameter, $\zeta_{t+1}\coloneqq\zeta_{1}\rho^{t}$ for
$t\geq 0$. We use softplus activations to ensure $\zeta_{0},\zeta_{1},\eta>0$
and the sigmoid function to make $0<\rho<1$.

It remains to specify the loss function used to learn the hyperparameters. One
immediate idea is to optimize the performance in reconstructing the true
tensor $\bm{\mathcal{X}}_{\star}$ by minimizing

$\displaystyle\mathcal{L}_{\text{SL}}(\bm{F})\coloneqq\frac{\left\|\bm{\mathcal{X}}_{\star}-\big{(}\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)}\big{)}\operatorname{\bm{\cdot}}\bm{\mathcal{G}}\right\|_{\operatorname{\mathsf{F}}}^{2}}{\left\lVert\bm{\mathcal{X}}_{\star}\right\rVert_{\operatorname{\mathsf{F}}}^{2}},$
(6)

which we refer to as supervised learning (SL) similar to [CLY21]. While this
is possible when trained on synthetic data—where we have access to the ground
truth—it may not be applicable in general in the absence of ground truth.

One might naturally wonder if (2) can be used to optimize the hyperparameters.
Unfortunately, with the decaying threshold schedule, (2) can be trivially
solved as we increase the number of iterations $T$ or choosing $\rho$ close to
$0$, since the threshold $\zeta_{t}$ will rapidly approach $0$, making
$\bm{\mathcal{S}}_{t}\approx\bm{\mathcal{Y}}-\bm{\mathcal{X}}_{t}$ by (5a).
Hence, the objective function (2) is no longer appropriate in learning the
hyperparameters. As such, we will instead use the following loss function:

$\displaystyle\mathcal{L}_{\text{SSL}}(\bm{F})\coloneqq\frac{1}{\left\lVert\bm{\mathcal{Y}}\right\rVert_{\operatorname{\mathsf{F}}}^{2}}\left\|\bm{\mathcal{Y}}-\big{(}\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)}\big{)}\operatorname{\bm{\cdot}}\bm{\mathcal{G}}\right\|_{1}$
(7)

to encourage sparsity of the corruption tensor [CCD+21, HKK+20]. We refer to
this loss function as the self-supervised learning (SSL) loss since it does
not require the ground truth and therefore, more amenable to real-data
scenarios.

###  3.2 From supervised learning to self-supervised learning

Depending on the loss function, our model can be used for both supervised or
self-supervised learning of the hyperparameters. For supervised learning,
suppose we have access to $\bm{\mathcal{X}}_{\star}$ for all tensor RPCA
problem instances in a training dataset, then the hyperparameters can be tuned
by optimizing (6).

A drawback of a purely supervised approach is the assumption of a “one size
fits all” set of hyperparameters that can be applied to other similar tensors.
In fact, tensor inputs into RPCA are exceptionally complex and in general,
cannot be adequately described in a few features such as their sparsity
levels, ranks, and condition numbers. Fortunately, the tuned hyperparameters
obtained from training serve as a warm start for further refinements. Similar
to [HKK+20], we first learn a fixed set of hyperparameters on the training
dataset (when available), which will be fine tuned for each tensor during test
time by minimizing (7) individually for each tensor, leading to a self-
supervised learning paradigm.

##  4 Experiments

We conduct synthetic and real-world data experiments to corroborate the
effectiveness of the proposed learned tensor RPCA approach. Code available at

https://github.com/hdong920/Tensor_RPCA_ScaledGD.

###  4.1 Synthetic data

We explore the effect of fine tuning with SSL for each combination of
$\alpha\in\\{0,0.1,\dots,0.9\\}$ and $r\in\\{10,20,\dots,70\\}$ with fixed
$n=100$. When a rank $(r,r,r)$ tensor $\bm{\mathcal{Y}}\in\mathbb{R}^{n\times
n\times n}$ needs to be sampled, we first randomly generate factor matrices
$\bm{U}_{\star}^{(1)},\bm{U}_{\star}^{(2)},\bm{U}_{\star}^{(3)}\in\mathbb{R}^{n\times
r}$ with orthonormal columns and core tensor $\bm{\mathcal{G}}_{\star}$ where
$[\bm{\mathcal{G}}_{\star}]_{i,i,i}=\kappa^{-(i-1)/(r-1)}$ for $\kappa=5$ and
$0$ elsewhere to construct
$\bm{\mathcal{X}}_{\star}=(\bm{U}_{\star}^{(1)},\bm{U}_{\star}^{(2)},\bm{U}_{\star}^{(3)})\operatorname{\bm{\cdot}}\bm{\mathcal{G}}_{\star}$.
In addition, $\bm{\mathcal{S}}_{\star}$ is a sparse tensor with
$\alpha$-fraction of its entries drawn from a uniform distribution in
$(-\theta,\theta)$, where
$\theta\coloneqq\frac{1}{n^{3}}\left\lVert\bm{\mathcal{X}}_{\star}\right\rVert_{1}$.
All experiments in this section were run for $T=100$ iterations of Algorithm
1.

First, we perform supervised learning by minimizing (6) similar to [CLY21] for
each combination of $\alpha$ and $r$. For each gradient step, we generate a
new sample (which is a tensor RPCA problem instance). We train for 1000
gradient steps with a decaying learning rate scheduler. Next, we perform SSL
for 500 gradient steps to fine tune the hyperparameters for 20 samples
generated for each combination of $\alpha$ and $r$ by minimizing (7). As the
baseline, we use Optuna [ASY+19], a Bayesian optimization-based hyperparameter
tuning package, to minimize (7) individually for 20 samples. We run Optuna for
500 iterations for each sample over a search space wide enough to contain the
learned hyperparameters.

(a) Baseline

(b) Supervised

(c) Fine tuning

Figure 2: The log mean relative recovery error,
$\frac{\left\lVert\bm{\mathcal{X}}_{\star}-\bm{\mathcal{X}}_{T}\right\rVert_{\operatorname{\mathsf{F}}}}{\left\lVert\bm{\mathcal{X}}_{\star}\right\rVert_{\operatorname{\mathsf{F}}}}$
for each combination of $\alpha$ and $r$ of Algorithm 1 using hyperparameters
found by the baseline, supervised learning, and supervised learning with self-
supervised fine tuning.

From Fig. 2, we observe that the baseline accurately recovers
$\bm{\mathcal{X}}_{\star}$ across a wide range of cases, though it can be
unstable in low rank and high sparsity scenarios. Supervised learning improves
the stability of these scenarios but performs slightly worse along the edge of
the phase transition. Although fine tuning adds more computation, we
simultaneously obtain more stable results than the baseline and gain a
significant improvement on recovering $\bm{\mathcal{X}}_{\star}$ in a broad
set of scenarios over supervised learning, especially along the edge of the
phase transition.

During self-supervised fine tuning, we track the percent change of the RPCA
hyperparameters’ fined tuned values relative to their final values from
supervised learning. Table 1 contains the quartiles of these percent changes
for samples that observed a more than a 90% reduction in recovery error after
fine tuning. From this, we see that RPCA hyperparameters can be quite
sensitive where tuning by grid search would require very fine sampling. Fine
tuning adapts to the sensitivity of RPCA hyperparameters which supervised
learning alone cannot address.

Quartile |  $\zeta_{0}$ |  $\zeta_{1}$ |  $\eta$ |  $\rho$
---|---|---|---|---
Q1 |  0.93% |  3.06% |  -14.03% |  3.36%
Q2 |  5.17% |  28.91% |  -3.38% |  6.20%
Q3 |  15.73% |  36.56% |  7.85% |  6.68%
Table 1: Quartiles of fine tuning parameter percent change for samples with
greater than $90\%$ reduction in recovery loss after fine tuning.

###  4.2 Background subtraction in video surveillance

Next, we apply learned tensor RPCA from scratch on the background subtraction
task using the Background Models Challenge real videos dataset [VCWL12], which
contains 9 videos of varying shapes. We use the first 6 to train two separate
models via supervised learning and SSL with $T=150$ for 15 epochs, which are
tested on the last 3 videos. Models assumed rank 1 along the time dimension
and full rank for all others. For SSL, it minimizes (7). However, for
supervised learning, as this dataset includes only ground truth binary
foreground masks (Boolean labels that indicate if a pixel is part of the
foreground), we cannot use (6) directly. Defining $\bm{\mathcal{M}}$ to be the
binary foreground mask of the same shape as $\bm{\mathcal{Y}}$, we know that
$\bm{\mathcal{S}}_{\star}$ and $\bm{\mathcal{M}}$ share the same support.
Letting $\odot$ be the Hadamard product, we use

$\displaystyle\mathcal{L}_{\text{SM}}(\bm{F})\coloneqq\frac{\left\|\left(\bm{\mathcal{Y}}-\big{(}\bm{U}^{(1)},\bm{U}^{(2)},\bm{U}^{(3)}\big{)}\operatorname{\bm{\cdot}}\bm{\mathcal{G}}\right)\odot(1-\bm{\mathcal{M}})\right\|_{\operatorname{\mathsf{F}}}^{2}}{\left\lVert\bm{\mathcal{Y}}\odot(1-\bm{\mathcal{M}})\right\rVert_{\operatorname{\mathsf{F}}}^{2}}$

for supervised learning since
$\bm{\mathcal{Y}}\odot(1-\bm{\mathcal{M}})=\bm{\mathcal{X}}_{\star}\odot(1-\bm{\mathcal{M}})$.
Because we only assume low rank along the time dimension, iterative updates of
Algorithm 1 are skipped for other modes to save computation. Fig. 3 display no
significant visual difference between results from the two models, suggesting
the benefit of SSL in alleviating the need of labeled training data for this
task. This is backed up by their similar masked average relative recovery test
errors, $\sqrt{\mathcal{L}_{\text{SM}}(\bm{F}_{T})}$. We observe that as the
assumed time dimension rank increased, the more fine tuning improved relative
recovery error.

(a) Supervised learning

(b) Self-supervised learning

Figure 3: Background subtraction of a video containing a moving train and
worker. From left to right of each half, the columns are the input video,
extracted sparse foreground, and low-rank background.

###  4.3 Electron backscatter diffraction microscopy

We compare our SSL-based tensor RPCA method and the baseline (tuned by Optuna
[ASY+19] using the loss (7)) on high-dimensional microscopy data from
[CSD+21]. To provide context, data collection involved iteratively performing
electron backscatter diffraction (EBSD) on a nickel-based superalloy. EBSD
collects an image where the orientation of the local crystal, represented as
Euler angles, is stored on each pixel. The superalloy was mechanically
polished after each EBSD image collection, removing approximately 1 micron of
material. This was repeated many times to obtain a sequence of slices
containing orientation measurements, with an in-plane pixel spacing of 1
micron. The result is a volume of orientation information that can be
considered as a high-dimensional tensor. See Fig. 4(a) for a visualization of
consecutive slices. Note this dataset has no inherent notion of ground truth
that we can perform supervised learning on for RPCA and contains only one
sample.

The motivation of applying RPCA here is that materials scientists usually
arduously parse through these these slices manually. Based on Fig. 4(a),
changes from slice to slice are very subtle. Being able to separate the large
slow moving sections from the sparse rapidly evolving parts could allow
materials scientists to quickly identify critical rare events with
implications on life-limiting behavior. Due to the dataset’s significant
amount of structure, we train our model for only 10 iterations of
hyperparameter updates on $(250\times 250\times 250\times 3)$ samples where we
identify unique parameters for each sample using $\bm{r}=(250,25,250,3)$ and
$T=10$, again skipping updates along full rank modes. From Fig. 4(b) and 4(c),
our method extracts sparse components that capture detailed changes across
slices and low-rank components that contain macro-scale patterns. We repeat
the experiment for 50 iterations of the baseline method, with the results
shown in Fig. 4(d) and 4(e). Although the baseline resulted in similar loss
values and captured similar quality low rank and sparse structures, it took
many more iterations (about 5 times more) to obtain these results.
Furthermore, we note that the per-iteration hyperparameter update runtimes of
the baseline and our model were very similar, so our model also ran about 5
times faster.

(a) Original slices

(b) Low rank component produced by SSL

(c) Sparse component produced by SSL

(d) Low rank component produced by Baseline

(e) Sparse component produced by Baseline

Figure 4: The tensor RPCA outputs on example EBSD slices.

##  5 Conclusions

We have illustrated and empirically shown two ways deep unfolded SSL can be
used to improve current approaches to tensor RPCA. First, a deep unfolded
version of Algorithm 1 with a new $\ell_{1}$-based loss function to train the
hyperparameters can perform equivalently with supervised models and quickly
provide satisfactory results in tasks where supervision is inapplicable.
Second, self-supervised fine tuning improves recovery by finding a unique set
of RPCA hyperparameters for each RPCA problem instance. Altogether, we have
simultaneously extended the applicability of tensor RPCA to data-starved
settings and enhanced its performance. Possible future adaptations include
extending our method to design similar models for tensor algorithms that are
centered around reconstruction like tensor completion. It would also be of
interest to replace our recurrent layers with other architectures designed for
sequential inputs, such as transformers.

## Acknowledgments

The work of H. Dong and Y. Chi is supported in part by the Air Force D3OM2S
Center of Excellence under FA8650-19-2-5209, by ONR under N00014-19-1-2404, by
the DARPA TRIAD program under Agreement No. HR00112190099, by NSF under
CCF-1901199 and ECCS-2126634, and by the Carnegie Mellon University
Manufacturing Futures Initiative, made possible by the Richard King Mellon
Foundation.

## References

  * [AJSN16] A. Anandkumar, P. Jain, Y. Shi, and U. N. Niranjan.  Tensor vs. matrix methods: Robust tensor decomposition under block sparse perturbations.  In Artificial Intelligence and Statistics, pages 268–276. PMLR, 2016.
  * [ASY+19] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2019.
  * [CCD+21] V. Charisopoulos, Y. Chen, D. Davis, M. Díaz, L. Ding, and D. Drusvyatskiy.  Low-rank matrix recovery with composite optimization: good conditioning and rapid convergence.  Foundations of Computational Mathematics, 21(6):1505–1593, 2021\.
  * [CLY21] H. Cai, J. Liu, and W. Yin.  Learned robust PCA: A scalable deep unfolding approach for high-dimensional outlier detection.  Advances in Neural Information Processing Systems, 34:16977–16989, 2021.
  * [CSD+21] M. G. Chapman, M. N. Shah, S. P. Donegan, J. M. Scott, P. A. Shade, D. Menasche, and M. D. Uchic.  Afrl additive manufacturing modeling series: challenge 4, 3d reconstruction of an in625 high-energy diffraction microscopy sample using multi-modal serial sectioning.  Integrating Materials and Manufacturing Innovation, 10(2):129–141, 2021.
  * [CZS+19] R. Cohen, Y. Zhang, O. Solomon, D. Toberman, L. Taieb, R. J. van Sloun, and Y. C. Eldar.  Deep convolutional robust PCA with application to ultrasound imaging.  In ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3212–3216. IEEE, 2019.
  * [DBBG19] D. Driggs, S. Becker, and J. Boyd-Graber.  Tensor robust principal component analysis: Better recovery with atomic norm regularization.  arXiv preprint arXiv:1901.10991, 2019.
  * [DTMC22] H. Dong, T. Tong, C. Ma, and Y. Chi.  Fast and provable tensor robust principal component analysis via scaled gradient descent.  arXiv preprint arXiv:2206.09109, 2022.
  * [FL18] S. Friedland and L.-H. Lim.  Nuclear norm of higher-order tensors.  Mathematics of Computation, 87(311):1255–1281, 2018.
  * [GL10] K. Gregor and Y. LeCun.  Learning fast approximations of sparse coding.  In Proceedings of the 27th international conference on international conference on machine learning, pages 399–406, 2010.
  * [HKK+20] C. Herrera, F. Krach, A. Kratsios, P. Ruyssen, and J. Teichmann.  Denise: Deep learning based robust PCA for positive semidefinite matrices.  stat, 1050:5, 2020.
  * [LFC+16] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan.  Tensor robust principal component analysis: Exact recovery of corrupted low-rank tensors via convex optimization.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5249–5257, 2016.
  * [LFC+19] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, and S. Yan.  Tensor robust principal component analysis with a new tensor nuclear norm.  IEEE transactions on pattern analysis and machine intelligence, 42(4):925–938, 2019.
  * [SCZ+19] O. Solomon, R. Cohen, Y. Zhang, Y. Yang, Q. He, J. Luo, R. J. van Sloun, and Y. C. Eldar.  Deep unfolded robust PCA with application to clutter suppression in ultrasound.  IEEE transactions on medical imaging, 39(4):1051–1063, 2019.
  * [TMC21a] T. Tong, C. Ma, and Y. Chi.  Accelerating ill-conditioned low-rank matrix estimation via scaled gradient descent.  Journal of Machine Learning Research, 22:1–63, 2021.
  * [TMC21b] T. Tong, C. Ma, and Y. Chi.  Low-rank matrix recovery with scaled subgradient methods: Fast and robust convergence without the condition number.  IEEE Transactions on Signal Processing, 69:2396–2409, 2021.
  * [TMC22] T. Tong, C. Ma, and Y. Chi.  Accelerating ill-conditioned robust low-rank tensor regression.  In ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 9072–9076. IEEE, 2022.
  * [TMPB+22] T. Tong, C. Ma, A. Prater-Bennette, E. Tripp, and Y. Chi.  Scaling and scalability: Provable nonconvex low-rank tensor estimation from incomplete measurements.  Journal of Machine Learning Research, 23(163):1–77, 2022.
  * [VCWL12] A. Vacavant, T. Chateau, A. Wilhelm, and L. Lequievre.  A benchmark dataset for outdoor foreground/background extraction.  In Asian Conference on Computer Vision, pages 291–300. Springer, 2012.
  * [VLJED21] H. Van Luong, B. Joukovsky, Y. C. Eldar, and N. Deligiannis.  A deep-unfolded reference-based RPCA network for video foreground-background separation.  In 2020 28th European Signal Processing Conference (EUSIPCO), pages 1432–1436. IEEE, 2021.
  * [YZ16] M. Yuan and C.-H. Zhang.  On tensor completion via nuclear norm minimization.  Foundations of Computational Mathematics, 16(4):1031–1068, 2016\.


# Divergence of Majorana-Phonon Scattering in Kitaev Quantum Spin Liquid

Haoxiang Li  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  A. Said  Advanced Photon Source,
Argonne National Laboratory, Argonne, Illinois 60439, USA  J.Q. Yan  Materials
Science and Technology Division, Oak Ridge National Laboratory, Oak Ridge,
Tennessee 37831, USA  D.M. Mandrus  Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Department of
Materials Science and Engineering, University of Tennessee at Knoxville,
Knoxville, Tennessee 37996, USA  H. N. Lee  Materials Science and Technology
Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  S.
Okamoto  Materials Science and Technology Division, Oak Ridge National
Laboratory, Oak Ridge, Tennessee 37831, USA  Gábor B. Halász<EMAIL_ADDRESS>Materials Science and Technology Division, Oak Ridge National Laboratory, Oak
Ridge, Tennessee 37831, USA  H. Miao<EMAIL_ADDRESS>Materials Science and
Technology Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee
37831, USA

(August 28, 2024)

###### Abstract

Magnetoelastic interaction couples spin and lattice degrees of freedom and
plays a key role in thermal transport properties of magnetic insulators. In
the Kitaev quantum spin liquid, the low energy excitations are charge neutral
Majorana fermions, which transform the magnetoelasctic interaction into
Majorana-phonon scattering. Motivated by anomalous thermal properties of the
Kitaev quantum spin liquid candidate $\alpha$-RuCl3, in this letter, we
combine meV resolution inelastic x-ray scattering and theoretical calculation
to examine the Majorana-phonon scattering. We analytically derive the
velocity-dependent Majorana-phonon scattering and find a divergence when the
acoustic phonons and the itinerant Majorana fermions have the same velocity.
Based on the experimentally determined acoustic phonon velocity in
$\alpha$-RuCl3, we estimate the range in the Kitaev interaction for which
divergent Majorana-phonon scattering can happen. Our result opens a new avenue
to uncover fractionalized quasiparticles in the Kitaev quantum spin liquid and
emphasizes the critical role of lattice excitations in $\alpha$-RuCl3.

Thermal transport is one of the key physical properties to uncover the nature
of emergent quantum states, such as unconventional superconductivity [1, 2, 3,
4, 5] and other exotic states that host anyonic excitations [6, 7, 8, 9, 10,
11, 12, 13, 14, 15, 16]. In the Kitaev honeycomb model [6],

$H=\sum_{\gamma,<i,j>}J^{\gamma}_{K}S^{\gamma}_{i}S^{\gamma}_{j}$ (1)

(where $J^{\gamma}_{K}$ ($\gamma=X,Y,Z$) is the bond-dependent coupling
parameter, and $<i,j>$ stands for nearest-neighbor pairs of spins at one of
the $X$, $Y$, or $Z$ bonds), the low-energy spin excitations fractionalize
into charge neutral Majorana fermions. Interestingly, theoretical studies
found that Majorana fermions can couple with collective lattice excitations,
known as phonons, through the magnetoelastic interaction and give rise to
unusual thermal transport behaviors [17, 18, 8].

$\alpha$-RuCl3, consisting of edge-shared Ru-Cl octahedra [Fig. 1(a)], is
widely believed as a promising Kitaev quantum spin liquid candidate [13, 14,
6, 19, 20, 21, 22, 23, 8, 24, 25, 26]. As shown in Fig. 1(b), the destructive
quantum-interference through the close-to-90∘ Ru–Cl–Ru bonds significantly
suppresses the Heisenberg magnetic exchange interaction, yielding a dominant
Ising-type interaction ($J_{K}$) perpendicular to the Ru–Cl–Ru plane [27].
Although non-Kitaev interactions are present in $\alpha$-RuCl3 [22, 24, 19,
28, 29], existence of Majorana fermions has been suggested by numerous
experimental studies [21, 22, 24, 30, 13, 25, 26, 14, 31, 15, 16, 32]. Most
interestingly, a half-integer thermal Hall conductivity in units of
$\pi^{2}k_{B}^{2}T/3h$ ($k_{B}$ is the Boltzmann constant) is observed [13,
14, 15], supporting non-Abelian anyons in $\alpha$-RuCl3 under in-plane
magnetic field. However, the quantized thermal Hall conductivity is found to
be sample-dependent [13, 14, 15, 16, 32], posting questions on the subtle
interplay between electronic excitations and lattice degrees of freedom [33].
In this letter, we combine meV-resolution inelastic x-ray scattering (IXS) and
analytical calculation to investigate the Majorana-phonon scattering process.
We find that the Majorana-phonon scattering rate diverges when the velocity of
acoustic phonons matches the velocity of Majorana fermions. Based on the IXS
determined acoustic phonon velocity, we identify the range in the Kitaev
interaction, $J_{K}$, for which divergent Majorana-phonon scattering can
happen in $\alpha$-RuCl3, and we find that the two velocities may be very
close for reasonable values of $J_{K}$. Our result indicates that structural
defects, such as stacking faults, that modify magnetic properties [22, 24],
may strongly affect the thermal transport behaviors in $\alpha$-RuCl3.

Millimeter-sized $\alpha$-RuCl3 crystals were grown by the sublimation of
$\alpha$-RuCl3 powder sealed in a quartz tube under vacuum [34]. Magnetic
order was confirmed to occur at 7 K by measuring magnetic properties and
specific heat [22]. The IXS experiments were conducted at beam line 30-ID-C
(HERIX) at the Advanced Photon Source (APS). The highly monochromatic X-ray
beam of incident energy $E_{i}=23.7$ keV ($\lambda=0.5226\leavevmode\nobreak\
\AA$) was focused on the sample with a beam cross section of $\sim 35\times
15$ $\mathrm{\mu m}^{2}$ (horizontal$\times$vertical). The total energy
resolution of the monochromatic X-ray beam and analyzer crystals was $\Delta
E\sim 1.5$ meV (full width at half maximum). The measurements were performed
in transmission geometry. Typical counting times were in the range of 30–120 s
per point in the energy scans at constant momentum transfer Q. H, K, L are
defined in the trigonal structure with $a=b=5.9639$ $\AA$, $c=17.17$ $\AA$ at
the room temperature.

Figure 1: Schematics of Majorana-phonon scattering. (a)-(b) Lattice structure
of $\alpha$-RuCl3. The red, blue and green lines in (a) represent the $X$, $Y$
and $Z$ bonds in the Kitaev model (see Eq. 1). (c)-(e) Schematics of the
scattering between phonons (green cone) and MFs (yellow cone) with
$v_{\mathrm{ph}}>v_{\mathrm{f}}$, $v_{\mathrm{ph}}=v_{\mathrm{f}}$ and
$v_{\mathrm{ph}}<v_{\mathrm{f}}$.

Before analyzing the IXS result, we analytically derive the velocity dependent
Majorana-phonon scattering in a gapless Kitaev spin liquid. Focusing on the
lowest temperatures, we consider acoustic phonons of velocity
$v_{\mathrm{ph}}$ and linearly dispersing Majorana fermions of velocity
$v_{\mathrm{f}}$; the two Dirac points of Majorana fermions at momenta
$\pm\mathbf{K}$ (which correspond to the $\mathrm{K}$ and
$\mathrm{K}^{\prime}$ points of the Brillouin zone) are equivalent to a single
Dirac point of complex fermions at momentum $\mathbf{K}$. The simplest
processes annihilating a phonon of momentum $\mathbf{q}$ then involve the
scattering of a fermion with relative momentum $\mathbf{k}$ (with respect to
$\mathbf{K}$) into another fermion with relative momentum
$\mathbf{k}+\mathbf{q}$. Note that processes creating or annihilating two
fermions are forbidden due to the nontrivial momentum $\mathbf{K}$ of the
fermion Dirac point ($2\mathbf{K}$ is not a reciprocal lattice vector). We
also remark that, while the phonons are generally three dimensional, only
those with a purely in-plane (two-dimensional) momentum $\mathbf{q}$ can
participate in the above scattering processes.

Importantly, due to the conservation of energy and momentum, the Majorana-
phonon scattering is strongly velocity dependent. For
$v_{\mathrm{ph}}>v_{\mathrm{f}}$ [see Fig. 1(c)], a negative-energy fermion
can be scattered into a positive-energy fermion, whereas for
$v_{\mathrm{ph}}<v_{\mathrm{f}}$ [see Fig. 1(e)], a positive-energy fermion
can scatter into a positive-energy fermion or a negative-energy fermion can
scatter into negative-energy fermion [17, 18]. Furthermore, as we show in the
Supplemental Material [35], the real and imaginary parts of the phonon self
energy (corresponding to the energy renormalization and the decay rate,
respectively) both diverge for $v_{\mathrm{ph}}\to v_{\mathrm{f}}$ [see Fig.
1(d)]. The divergent components of these two quantities at momentum
$\mathbf{q}$ and temperature $T=1/\beta$ are found to be

$\mathrm{Re}\Sigma_{\mathbf{q}}\propto\bigg{\\{}\begin{array}[]{c}+|\mathbf{q}|^{2}\left(v_{\mathrm{ph}}-v_{\mathrm{f}}\right)^{-1/2}\mathcal{I}^{(2)}(\beta
v_{\mathrm{f}}|\mathbf{q}|)\quad(v_{\mathrm{ph}}>v_{\mathrm{f}}),\\\
-|\mathbf{q}|^{2}\left(v_{\mathrm{f}}-v_{\mathrm{ph}}\right)^{-1/2}\mathcal{I}^{(1)}(\beta
v_{\mathrm{f}}|\mathbf{q}|)\quad(v_{\mathrm{ph}}<v_{\mathrm{f}}),\end{array}$
(2)
$\mathrm{Im}\Sigma_{\mathbf{q}}\propto\bigg{\\{}\begin{array}[]{c}|\mathbf{q}|^{2}\left(v_{\mathrm{ph}}-v_{\mathrm{f}}\right)^{-1/2}\mathcal{I}^{(1)}(\beta
v_{\mathrm{f}}|\mathbf{q}|)\quad(v_{\mathrm{ph}}>v_{\mathrm{f}}),\\\
|\mathbf{q}|^{2}\left(v_{\mathrm{f}}-v_{\mathrm{ph}}\right)^{-1/2}\mathcal{I}^{(2)}(\beta
v_{\mathrm{f}}|\mathbf{q}|)\quad(v_{\mathrm{ph}}<v_{\mathrm{f}}),\end{array}$
(3)

where $\mathcal{I}^{(1,2)}(x)$ with $x=\beta v_{\mathrm{f}}|\mathbf{q}|$ are
positive dimensionless functions with asymptotic forms

$\mathcal{I}^{(1)}(x)\sim\bigg{\\{}\begin{array}[]{c}x\quad(x\ll 1),\\\
1\quad(x\gg
1),\end{array}\quad\mathcal{I}^{(2)}(x)\sim\bigg{\\{}\begin{array}[]{c}x^{-1}\quad\,\,\,\,(x\ll
1),\\\ x^{-3/2}\quad(x\gg 1).\end{array}$ (4)

Hence, the real and imaginary parts of the phonon self energy have completely
different scaling behaviors as a function of $\beta
v_{\mathrm{f}}|\mathbf{q}|$. In particular, at the smallest phonon momenta
($\beta v_{\mathrm{f}}|\mathbf{q}|\ll 1$), the divergent contribution is
expected to be most observable in the real part for
$v_{\mathrm{ph}}>v_{\mathrm{f}}$ and in the imaginary part for
$v_{\mathrm{ph}}<v_{\mathrm{f}}$. We also note that the plus (minus) sign in
Eq. 2 corresponds to hardening (softening) of the phonon energy.

Figure 2: Acoustic phonon excitations at small q. (a) Schematics of the IXS
experimental setup and trajectory of the measurement in Brillouin zone.
(b)-(d) Three branches of acoustic phonon excitations measured by IXS, where
LA is measured from the Q=(6, -3, 0)–(6.2, -3.1, 0), TA is measured from Q=(6,
-3, 0)–(6, -2.75, 0), TAO is measured from Q=(6, -3, 0) – (6, -3, 0.266) The
plot shows the Bose-factor corrected IXS intensity. (e)-(g) Raw IXS spectra
(solid circles) of the three phonon branches and the fits to the data (dotted
lines). The error bars represent one standard deviation assuming Poisson
counting statistics.

We now determine the in-plane acoustic phonon velocities using IXS. IXS
directly probes the phonon dynamical structure factor,
$S(\mathbf{Q=q+G},\omega)=\chi\prime\prime(\mathbf{Q},\omega)/(1-e^{-\omega/k_{B}T})$,
where $\mathbf{Q}$ and $\mathbf{q}$ are total and reduced momentum transfer.
$\mathbf{G}$ is the reciprocal lattice vector and $\omega$ is the phonon
energy. $\chi\prime\prime(\mathbf{Q},\omega)$ is the the imaginary part of the
dynamical phonon susceptibility that can be obtained from
$S(\mathbf{Q=q+G},\omega)$. For acoustic modes, all atoms in the unit cell
move in the same direction as the polarization vector,
$\mathbf{e}(\mathbf{q})$, and give rise to
$S(\mathbf{Q},\omega)\propto|\mathbf{Q}\cdot\mathbf{e}(\mathbf{q})|^{2}$.
Figure 2(a) shows the experimental geometry in the trigonal structure at room
temperature. To selectively enhance the in-plane longitudinal (LA) and
transverse (TA) acoustic phonons, we perform IXS measurement along (6, -3,
0)$\longrightarrow$(7, -3.5, 0) and (6, -3, 0)$\longrightarrow$(6, -2, 0) in
reciprocal lattice unit (r.l.u.), respectively. We also performed measurement
along (6, -3, 0)$\longrightarrow$(6, -3, 1) to evaluate the out-of-plane
transverse acoustic phonon (TAo).

Figure 2 shows the image plots of Bose factor corrected dynamical structure
factor of in-plane LA [Fig. 2(b)], in-plane TA [Fig. 2(c)] and out-of-plane
TAo [Fig. 2(d)]. Near the Brillouin zone center, the acoustic phonon velocity
is defined as $v_{\mathrm{ph}}=\frac{\partial\omega}{\partial
q}|_{q\rightarrow 0}\sim\frac{\omega}{q}$. In our measurement, the LA mode has
the largest slope at the Brillouin zone center and hence the largest
$v_{\mathrm{ph}}$. The velocity of TAo is the smallest at the Brillouin zone
center due to the van der Waals bonding between quasi-two dimensional layers.
Indeed, at $L=0.5$ r.l.u., the TAo mode, corresponding to a shear displacement
between adjacent $\alpha$-RuCl3-layers, has an energy of 3 meV at 300 K (see
supplemental materials [35]), consistent with large stacking fault in
$\alpha$-RuCl3. To precisely extract $\omega(\mathbf{q})$, we fit the raw IXS
spectra shown in Figs. 2(e)-(g). Phonon modes are represented by standard
damped harmonic oscillator functions (DHOFs) [36]:

$\mathrm{DHOF}=\sum_{i}I_{i}[\frac{\Gamma_{i}}{(\omega-\omega_{\mathbf{Q},i})^{2}+\Gamma_{i}^{2}}-\frac{\Gamma_{i}}{(\omega+\omega_{\mathbf{Q},i})^{2}+\Gamma_{i}^{2}}]$
(5)

where $I_{i}$, $\omega_{\mathbf{Q},i}$, and $2\Gamma_{i}$ are the intensity,
energy, and full width at half maximum (FWHM) of phonon peak $i$. To fit to
the IXS spectra, the DHOFs are weighted by the Bose factor and convoluted with
the 1.5 meV energy resolution. Details of the spectral peak fitting are given
in supplemental materials [35]. The extracted phonon peak dispersions are
shown in Fig. 3(a). The red, orange and blue symbols are corresponding to LA,
TA and TAo, respectively. The dashed lines are sinusoidal fittings of the
extracted $\omega(\mathbf{q})$.

As shown in Eqs. 2 and 3, the Majorana-phonon scattering diverges when
$v_{\mathrm{ph}}=v_{\mathrm{f}}$. In the isotropic limit of the pure Kitaev
model with $J_{K}^{X}=J_{K}^{Y}=J_{K}^{Z}=J_{K}$, the Majorana fermion
dispersion near the $K$ point is
$\omega(\mathbf{k})=\frac{\sqrt{3}}{4}J_{K}a|\mathbf{k}|$, where the in-plane
lattice constant is $a=5.9639$ Å. The experimentally determined
$v_{\mathrm{ph}}$ can thus be used to identify the range in $J_{K}$ for which
divergent Majorana-phonon scattering can happen. Using the largest acoustic
phonon velocity, $v^{\mathrm{LA}}_{\mathrm{ph}}=23$ meV$\cdot$Å, we obtain
$J^{max}_{K}\sim 9$ meV.

Previous experimental and theoretical studies derived a wide range of
$|J_{K}|$ from 1 to 25 meV, although most studies support $5<|J_{K}|<13$ meV.
In the inelastic neutron scattering study, $J_{K}\sim$ 7 meV was extracted by
fitting the spin-wave with a model Hamiltonian [22], which corresponds to
$v_{\mathrm{f}}\sim 18$ meV$\cdot$Å. This value is very close to the phonon
velocity of the TA mode with $v^{\mathrm{TA}}_{\mathrm{ph}}=17$ meV$\cdot$Å,
suggesting a close to divergent scattering. Interestingly, a recent IXS study
found a 15% phonon softening in the TA mode at $q\leq$0.1 r.l.u. [31, 35].
This softening effect is consistent with the Majorana-phonon scattering
induced phonon anomaly described by the real part of the phonon self energy in
Eq. 2. It is important to emphasize that $J_{K}\sim$ 7 meV may be incompatible
with the high fractionalization temperature, $T^{*}\sim$ 80 K, observed in
thermal transport measurements [13]. Numerical study of a pure Kitaev model
shows that $T^{*}\sim 0.5J_{K}$ [8]. For $J_{K}\sim$ 7 meV, the estimated
$T^{*}\sim 40$ K is only half of the experimental value.

According to Eq. 3, the phonon softening will also be accompanied by phonon
broadening at low temperature. It is interesting to compare the phonon
broadening effect in the electron-phonon coupled charge density wave (CDW)
systems, such as the cuprate high-Tc superconductors [37, 38] and transition
metal dichalcogenides [39]. The phonon broadening is usually measured by the
damping ratio, $\Gamma/\omega_{0}$, where $\Gamma$ is the full width at half
maximum of the phonon peak, $\omega_{0}$ is the bare phonon energy. For a 15 %
phonon softening at $\omega_{0}=2$ meV, as those observed in $\alpha$-RuCl3 at
10 K [31], the typical damping ratio in the CDW systems are about 0.1 to 0.3
[37, 39, 38], corresponding to $\Gamma\sim$ 0.2 to 0.6 meV. This broadening
effect is beyond the resolving power of the current study (energy resolution
$\Delta E\sim$ 1.5 meV). A similar case is the CDW in ZrTe3 [40], where the
phonon broadening effect at the low energy phonon peak $\omega_{0}\sim 2$ meV
was not resolved. A sub-meV resolution measurement, such as inelastic neutron
scattering, is required to uncover the Majorana-phonon scattering induced
phonon broadening.

Figure 3: Extracted phonon dispersions and velocities. (a) Phonon peak
dispersions extracted from the measured $S(\mathbf{Q},\omega)$. The x-axis is
in unit of inverse angstrom ($1/\mathrm{\AA}$), thus we can directly compare
the dispersion of three acoustic modes. The colored dash line is the
sinusoidal fit to the data. (b) The extracted phonon velocities at the
$\Gamma$ point based on the dispersions shown in panel (a). Error bars in
panel (a) represent the $2\sigma$ returned from spectral peak fittings.

Finally we discuss the implications of our result to the sample dependent
thermal transport measurement in $\alpha$-RuCl3 [13, 14, 15, 16, 32]. Due to
the van der Waals structure, the magnetoelastic coupling is expected to be
sensitive to material defects, such as stacking faults. Indeed, experimental
studies have found that the antiferromagnetic transition temperature can be
changed up to 50% depending on the level of stacking faults. Taking the
experimentally estimated $|J_{K}|$ between 5 and 13 meV, our results indicate
that the Majorana-phonon scattering in $\alpha$-RuCl3 may be close to the
divergent condition. In this case, the subtle change in structural and/or
magnetic interactions can amplify the Majorana-phonon scattering, giving rise
to a large impact on their thermal transport properties.

In summary, we investigated the velocity-dependent Majorana-phonon scattering
in the Kitaev quantum spin liquid and found that it shows divergent behavior
when the phonon and Majorana velocities match. Our result helps to establish
the phonon dynamical structure factor as a promising avenue to uncover
fractionalized excitations and sheds light on the novel thermal transport
properties of the Kiteav spin liquid candidate $\alpha$-RuCl3.

Acknowledgement: We thank T. Berlijn, H. Ding, J. K. Keum, G. Kotliar, S.
Nagler, N. Perkins, and A. Tennant for stimulating discussions. This research
at Oak Ridge National Laboratory (ORNL) was sponsored by the U.S. Department
of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and
Engineering Division (theory, IXS experiment and material synthesis). This
research used resources of the Advanced Photon Source, a U.S. Department of
Energy (DOE) Office of Science User Facility, operated for the DOE Office of
Science by Argonne National Laboratory under Contract No. DE-AC02-06CH11357.
Extraordinary facility operations were supported, in part, by the DOE Office
of Science through the National Virtual Biotechnology Laboratory, a consortium
of DOE national laboratories focused on the response to COVID-19, with funding
provided by the Coronavirus CARES Act. G.B.H. and S.O. were supported by the
U.S. Department of Energy, Office of Science, National Quantum Information
Science Research Centers, Quantum Science Center.

##  Appendix A Velocity dependence of the Majorana-phonon scattering

We consider the low-temperature scattering between acoustic phonons of
velocity $v_{\mathrm{ph}}$ and Majorana fermions of velocity $v_{\mathrm{f}}$
in a gapless Kitaev spin liquid. The two Dirac points of Majorana fermions at
momenta $\pm\mathbf{K}$ (which correspond to the $\mathrm{K}$ and
$\mathrm{K}^{\prime}$ points of the Brillouin zone) are equivalent to a single
Dirac point of complex fermions at momentum $\mathbf{K}$. The simplest
processes annihilating a phonon of momentum $\mathbf{q}$ then involve the
scattering of a fermion with relative momentum $\mathbf{k}$ (with respect to
$\mathbf{K}$) into another fermion with relative momentum
$\mathbf{k}+\mathbf{q}$. Note that the processes creating or annihilating two
fermions are forbidden due to the nontrivial momentum $\mathbf{K}$ of the
fermion Dirac point ($2\mathbf{K}$ is not a reciprocal lattice vector). Since
the Majorana-phonon coupling is proportional to the strain
$\epsilon_{\mathbf{q}}\sim|\mathbf{q}|u_{\mathbf{q}}$ with the lattice
displacement
$u_{\mathbf{q}}\sim(v_{\mathrm{ph}}|\mathbf{q}|)^{-1/2}(a_{\mathbf{q}}^{\phantom{{\dagger}}}+a_{-\mathbf{q}}^{{\dagger}})$
in terms of the bosonic phonon operators
$a_{\mathbf{q}}^{\phantom{{\dagger}}}$ and $a_{\mathbf{q}}^{{\dagger}}$, the
Majorana-phonon matrix element is proportional to
$(|\mathbf{q}|/v_{\mathrm{ph}})^{1/2}$. Hence, the imaginary part of the
phonon self energy (i.e., the phonon decay rate) at momentum $\mathbf{q}$,
energy $\omega$, and temperature $T$ is given by

$\displaystyle\mathrm{Im}\Sigma(\mathbf{q},\omega)$ $\displaystyle\propto$
$\displaystyle\frac{|\mathbf{q}|}{v_{\mathrm{ph}}}\int
d^{2}\mathbf{k}\,\bigg{\\{}\left[\tanh\left(\frac{v_{\mathrm{f}}|\mathbf{k}|}{2T}\right)+\tanh\left(\frac{v_{\mathrm{f}}|\mathbf{k}+\mathbf{q}|}{2T}\right)\right]\delta\big{[}\omega-
v_{\mathrm{f}}\left(|\mathbf{k}|+|\mathbf{k}+\mathbf{q}|\right)\big{]}$ (6)
$\displaystyle+2\left[\tanh\left(\frac{v_{\mathrm{f}}|\mathbf{k}+\mathbf{q}|}{2T}\right)-\tanh\left(\frac{v_{\mathrm{f}}|\mathbf{k}|}{2T}\right)\right]\delta\big{[}\omega-
v_{\mathrm{f}}\left(|\mathbf{k}+\mathbf{q}|-|\mathbf{k}|\right)\big{]}\bigg{\\}},$

where the first term describes the scattering of negative-energy fermions into
positive-energy fermions, while the second term accounts for the scattering of
positive-energy fermions into positive-energy fermions and the scattering of
negative-energy fermions into negative-energy fermions. We remark that, while
the phonons may generally be three dimensional, only those with a purely in-
plane (two-dimensional) momentum $\mathbf{q}$ can participate in the above
processes.

Focusing on the first term of Eq. (6), the delta function $\delta[\omega-
v_{\mathrm{f}}(|\mathbf{k}|+|\mathbf{k}+\mathbf{q}|)]$ necessarily vanishes
for $\omega<v_{\mathrm{f}}|\mathbf{q}|$ due to
$|\mathbf{k}|+|\mathbf{k}+\mathbf{q}|\geq|\mathbf{q}|$. In contrast, for
$\omega\geq v_{\mathrm{f}}|\mathbf{q}|$, the delta function is nonzero along
an ellipse defined by
$|\mathbf{k}|+|\mathbf{k}+\mathbf{q}|=\alpha|\mathbf{q}|$, where
$\alpha=\omega/(v_{\mathrm{f}}|\mathbf{q}|)\geq 1$ is the inverse
eccentricity. Taking $\mathbf{q}=(q,0)$ without loss of generality, and
introducing the elliptical coordinates
$\mathbf{k}=q(\cosh\mu\cos\nu-1,\sinh\mu\sin\nu)/2$, this ellipse corresponds
to a contour satisfying $\cosh\mu=\alpha$. Given that
$|\mathbf{k}|=q(\cosh\mu-\cos\nu)/2$ and
$|\mathbf{k}+\mathbf{q}|=q(\cosh\mu+\cos\nu)/2$, while
$d^{2}\mathbf{k}=q^{2}d\mu\,d\nu\,(\cosh^{2}\mu-\cos^{2}\nu)/4$ and

$\delta\big{[}\omega-
v_{\mathrm{f}}\left(|\mathbf{k}|+|\mathbf{k}+\mathbf{q}|\right)\big{]}=\delta\left[qv_{\mathrm{f}}\left(\alpha-\cosh\mu\right)\right]=\frac{\delta\left[\alpha-\cosh\mu\right]}{qv_{\mathrm{f}}}=\frac{\delta\left[\mu-\mathrm{arccosh}\,\alpha\right]}{qv_{\mathrm{f}}|\sinh\mu|}=\frac{\delta\left[\mu-\mathrm{arccosh}\,\alpha\right]}{qv_{\mathrm{f}}\sqrt{\alpha^{2}-1}},$
(7)

the phonon decay rate from the first term of Eq. (6) then becomes

$\displaystyle\mathrm{Im}\Sigma^{(1)}(\mathbf{q},\omega)$
$\displaystyle\propto$
$\displaystyle\frac{q^{2}}{4v_{\mathrm{f}}v_{\mathrm{ph}}\sqrt{\alpha^{2}-1}}\int_{0}^{2\pi}d\nu\left(\alpha^{2}-\cos^{2}\nu\right)\sum_{\xi=\pm}\tanh\left[\frac{qv_{\mathrm{f}}(\alpha+\xi\cos\nu)}{4T}\right]$
(8) $\displaystyle=$
$\displaystyle\frac{\sqrt{2}q^{3}}{v_{\mathrm{f}}}\left(\omega^{2}-q^{2}v_{\mathrm{f}}^{2}\right)^{-1/2}\mathcal{I}^{(1)}\left[\frac{\omega}{qv_{\mathrm{f}}},\frac{qv_{\mathrm{f}}}{T}\right]\quad(\omega>qv_{\mathrm{f}}),$

where the dimensionless function $\mathcal{I}^{(1)}[\alpha,\gamma]$ with
$\alpha=\omega/(qv_{\mathrm{f}})$ and $\gamma=qv_{\mathrm{f}}/T$ is given by
the integral

$\mathcal{I}^{(1)}[\alpha,\gamma]=\frac{1}{4\sqrt{2}\alpha}\int_{0}^{2\pi}d\nu\left(\alpha^{2}-\cos^{2}\nu\right)\sum_{\xi=\pm}\tanh\left[\frac{\gamma(\alpha+\xi\cos\nu)}{4}\right].$
(9)

Since the function $\mathcal{I}^{(1)}[\alpha,\gamma]$ is finite and well
behaved for all $\alpha\geq 1$ and $\gamma>0$, the phonon decay rate
$\mathrm{Im}\Sigma^{(1)}(\mathbf{q},\omega)$ from the first term of Eq. (6) is
divergent in the limit of $\alpha=\omega/(qv_{\mathrm{f}})\to 1$ due to the
factor $(\omega^{2}-q^{2}v_{\mathrm{f}}^{2})^{-1/2}$ in Eq. (8).

Focusing on the second term of Eq. (6), the delta function $\delta[\omega-
v_{\mathrm{f}}(|\mathbf{k}+\mathbf{q}|-|\mathbf{k}|)]$ necessarily vanishes
for $\omega>v_{\mathrm{f}}|\mathbf{q}|$ due to
$|\mathbf{k}+\mathbf{q}|-|\mathbf{k}|\leq|\mathbf{q}|$. In contrast, for
$\omega\leq v_{\mathrm{f}}|\mathbf{q}|$, the delta function is nonzero along a
hyperbola defined by
$|\mathbf{k}+\mathbf{q}|-|\mathbf{k}|=\alpha|\mathbf{q}|$, where
$\alpha=\omega/(v_{\mathrm{f}}|\mathbf{q}|)\leq 1$ is the inverse
eccentricity. Taking $\mathbf{q}=(q,0)$ without loss of generality, and
introducing the elliptical coordinates
$\mathbf{k}=q(\cosh\mu\cos\nu-1,\sinh\mu\sin\nu)/2$, this hyperbola
corresponds to a contour satisfying $\cos\nu=\alpha$. Given that
$|\mathbf{k}|=q(\cosh\mu-\cos\nu)/2$ and
$|\mathbf{k}+\mathbf{q}|=q(\cosh\mu+\cos\nu)/2$, while
$d^{2}\mathbf{k}=q^{2}d\mu\,d\nu\,(\cosh^{2}\mu-\cos^{2}\nu)/4$ and

$\delta\big{[}\omega-
v_{\mathrm{f}}\left(|\mathbf{k}+\mathbf{q}|-|\mathbf{k}|\right)\big{]}=\delta\left[qv_{\mathrm{f}}\left(\alpha-\cos\nu\right)\right]=\frac{\delta\left[\alpha-\cos\nu\right]}{qv_{\mathrm{f}}}=\frac{\delta\left[\nu-\mathrm{arccos}\,\alpha\right]}{qv_{\mathrm{f}}|\sin\nu|}=\frac{\delta\left[\nu-\mathrm{arccos}\,\alpha\right]}{qv_{\mathrm{f}}\sqrt{1-\alpha^{2}}},$
(10)

the phonon decay rate from the second term of Eq. (6) then becomes

$\displaystyle\mathrm{Im}\Sigma^{(2)}(\mathbf{q},\omega)$
$\displaystyle\propto$
$\displaystyle\frac{q^{2}}{2v_{\mathrm{f}}v_{\mathrm{ph}}\sqrt{1-\alpha^{2}}}\int_{-\infty}^{\infty}d\mu\left(\cosh^{2}\mu-\alpha^{2}\right)\sum_{\xi=\pm}\xi\tanh\left[\frac{qv_{\mathrm{f}}(\cosh\mu+\xi\alpha)}{4T}\right]$
(11) $\displaystyle=$
$\displaystyle\frac{\sqrt{2}q^{3}}{v_{\mathrm{f}}}\left(q^{2}v_{\mathrm{f}}^{2}-\omega^{2}\right)^{-1/2}\mathcal{I}^{(2)}\left[\frac{\omega}{qv_{\mathrm{f}}},\frac{qv_{\mathrm{f}}}{T}\right]\quad(\omega<qv_{\mathrm{f}}),$

where the dimensionless function $\mathcal{I}^{(2)}[\alpha,\gamma]$ with
$\alpha=\omega/(qv_{\mathrm{f}})$ and $\gamma=qv_{\mathrm{f}}/T$ is given by
the integral

$\mathcal{I}^{(2)}[\alpha,\gamma]=\frac{1}{2\sqrt{2}\alpha}\int_{-\infty}^{\infty}d\mu\left(\cosh^{2}\mu-\alpha^{2}\right)\sum_{\xi=\pm}\xi\tanh\left[\frac{\gamma(\cosh\mu+\xi\alpha)}{4}\right].$
(12)

Since the function $\mathcal{I}^{(2)}[\alpha,\gamma]$ is finite and well
behaved for all $0<\alpha\leq 1$ and $\gamma>0$, the phonon decay rate
$\mathrm{Im}\Sigma^{(2)}(\mathbf{q},\omega)$ from the second term of Eq. (6)
is divergent in the limit of $\alpha=\omega/(qv_{\mathrm{f}})\to 1$ due to the
factor $(q^{2}v_{\mathrm{f}}^{2}-\omega^{2})^{-1/2}$ in Eq. (11).

Figure S4: Dimensionless scaling functions $\mathcal{I}^{(1)}[\gamma]$ and
$\mathcal{I}^{(2)}[\gamma]$ (with $\gamma=qv_{\mathrm{f}}/T\approx
qv_{\mathrm{ph}}/T$) that govern the divergent components of the phonon decay
rate and the phonon energy renormalization in the limit of matching phonon and
fermion velocities ($v_{\mathrm{ph}}\to v_{\mathrm{f}}$).

In the limit of $\omega\to qv_{\mathrm{f}}$, the divergent component of the
net phonon decay rate,
$\mathrm{Im}\Sigma(\mathbf{q},\omega)=\mathrm{Im}\Sigma^{(1)}(\mathbf{q},\omega)+\mathrm{Im}\Sigma^{(2)}(\mathbf{q},\omega)$,
can be extracted by setting $\alpha=\omega/(qv_{\mathrm{f}})=1$ in the first
argument of each dimensionless function $\mathcal{I}^{(1,2)}[\alpha,\gamma]$.
Introducing the relative phonon energy $\tilde{\omega}=\omega-
qv_{\mathrm{f}}$, this divergent component takes the form

$\mathrm{Im}\tilde{\Sigma}(\mathbf{q},\tilde{\omega})\propto\sqrt{\frac{q^{5}}{v_{\mathrm{f}}^{3}|\tilde{\omega}|}}\left\\{\theta\left(\tilde{\omega}\right)\mathcal{I}^{(1)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}+\theta\left(-\tilde{\omega}\right)\mathcal{I}^{(2)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\right\\},$
(13)

where $\theta(\tilde{\omega})$ is the Heaviside step function, while the
dimensionless functions $\mathcal{I}^{(1)}[\gamma]$ and
$\mathcal{I}^{(2)}[\gamma]$ are given by

$\displaystyle\mathcal{I}^{(1)}[\gamma]\equiv\mathcal{I}^{(1)}[1,\gamma]=\frac{1}{4\sqrt{2}}\int_{0}^{2\pi}d\nu\sin^{2}\nu\sum_{\xi=\pm}\tanh\left[\frac{\gamma(1+\xi\cos\nu)}{4}\right],$
$\displaystyle\mathcal{I}^{(2)}[\gamma]\equiv\mathcal{I}^{(2)}[1,\gamma]=\frac{1}{2\sqrt{2}}\int_{-\infty}^{\infty}d\mu\sinh^{2}\mu\sum_{\xi=\pm}\xi\tanh\left[\frac{\gamma(\cosh\mu+\xi)}{4}\right].$
(14)

The precise forms of these functions are plotted in Fig. S4, while their
asymptotic behaviors in the limits of small and large $\gamma$ are

$\mathcal{I}^{(1)}[\gamma]\sim\Big{\\{}\begin{array}[]{c}\gamma\quad(\gamma\ll
1),\\\ 1\quad(\gamma\gg
1),\end{array}\qquad\mathcal{I}^{(2)}[\gamma]\sim\Big{\\{}\begin{array}[]{c}\gamma^{-1}\quad\,\,\,\,(\gamma\ll
1),\\\ \gamma^{-3/2}\quad(\gamma\gg 1).\end{array}$ (15)

Using the Kramers-Kronig relations, the divergent component of the phonon
energy renormalization, corresponding to the real part of the phonon self
energy, can then be obtained as

$\displaystyle\mathrm{Re}\tilde{\Sigma}(\mathbf{q},\tilde{\omega})$
$\displaystyle\propto$
$\displaystyle\frac{1}{\pi}\,\mathbb{P}\int_{-\infty}^{\infty}d\tilde{\omega}^{\prime}\,\frac{\mathrm{Im}\tilde{\Sigma}(\mathbf{q},\tilde{\omega}^{\prime})}{\tilde{\omega}-\tilde{\omega}^{\prime}}=\frac{1}{\pi}\sqrt{\frac{q^{5}}{v_{\mathrm{f}}^{3}}}\Bigg{\\{}\mathcal{I}^{(1)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\,\mathbb{P}\int_{0}^{\infty}\frac{d\tilde{\omega}^{\prime}}{\sqrt{\tilde{\omega}^{\prime}}(\tilde{\omega}-\tilde{\omega}^{\prime})}+\mathcal{I}^{(2)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\,\mathbb{P}\int_{0}^{\infty}\frac{d\tilde{\omega}^{\prime\prime}}{\sqrt{\tilde{\omega}^{\prime\prime}}(\tilde{\omega}+\tilde{\omega}^{\prime\prime})}\Bigg{\\}}$
(16) $\displaystyle=$
$\displaystyle\sqrt{\frac{q^{5}}{v_{\mathrm{f}}^{3}|\tilde{\omega}|}}\left\\{\theta\left(\tilde{\omega}\right)\mathcal{I}^{(2)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}-\theta\left(-\tilde{\omega}\right)\mathcal{I}^{(1)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\right\\}.$

Finally, by setting $\omega=qv_{\mathrm{ph}}$ and, hence,
$\tilde{\omega}=q(v_{\mathrm{ph}}-v_{\mathrm{f}})$, the physical decay rate
and energy renormalization are found to be

$\displaystyle\mathrm{Im}\tilde{\Sigma}(\mathbf{q})$ $\displaystyle\propto$
$\displaystyle\frac{q^{2}}{\sqrt{v_{\mathrm{f}}^{3}|v_{\mathrm{ph}}-v_{\mathrm{f}}|}}\left\\{\theta\left(v_{\mathrm{ph}}-v_{\mathrm{f}}\right)\mathcal{I}^{(1)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}+\theta\left(v_{\mathrm{f}}-v_{\mathrm{ph}}\right)\mathcal{I}^{(2)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\right\\},$
$\displaystyle\mathrm{Re}\tilde{\Sigma}(\mathbf{q})$ $\displaystyle\propto$
$\displaystyle\frac{q^{2}}{\sqrt{v_{\mathrm{f}}^{3}|v_{\mathrm{ph}}-v_{\mathrm{f}}|}}\left\\{\theta\left(v_{\mathrm{ph}}-v_{\mathrm{f}}\right)\mathcal{I}^{(2)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}-\theta\left(v_{\mathrm{f}}-v_{\mathrm{ph}}\right)\mathcal{I}^{(1)}\bigg{[}\frac{qv_{\mathrm{f}}}{T}\bigg{]}\right\\}.$
(17)

Both the phonon decay rate and the phonon energy renormalization diverge in
the limit of matching phonon and fermion velocities ($v_{\mathrm{ph}}\to
v_{\mathrm{f}}$). Moreover, the divergent contributions have specific scaling
behaviors as a function of $\gamma=qv_{\mathrm{f}}/T\approx
qv_{\mathrm{ph}}/T$ that can be observed by varying the phonon momentum
$\mathbf{q}$ and/or the temperature $T$. Interestingly, however, there are two
completely different scaling behaviors (see Fig. S4), and the one
characterizing the decay rate (energy renormalization) for
$v_{\mathrm{ph}}>v_{\mathrm{f}}$ characterizes the energy renormalization
(decay rate) for $v_{\mathrm{ph}}<v_{\mathrm{f}}$. Also, the divergent
component of the energy renormalization is positive for
$v_{\mathrm{ph}}>v_{\mathrm{f}}$ (corresponding to phonon hardening) but
negative for $v_{\mathrm{ph}}<v_{\mathrm{f}}$ (corresponding to phonon
softening).

##  Appendix B Inelastic X-ray scattering and spectral peak fitting

Inelastic X-ray (IXS) scattering is powerful tool for probing phonon
excitations with meV energy resolution [41]. In this work, the experimental
setup provides an energy resolution of 1.3 meV and a small X-ray beam size of
$\mathrm{35}\times\mathrm{15}$ $\mathrm{\mu m}^{2}$ . The energy resolution is
calibrated by fitting the elastic peak (Fig. S5) to a pseudo-voigt function:

$R(\omega)=(1-\alpha)\frac{I}{\sqrt{2\pi}\sigma}e^{-\frac{\omega}{2\sigma^{2}}}+\alpha\frac{I}{\pi}\frac{\Gamma}{\omega^{2}+\Gamma^{2}}$
(18)

where the energy resolution is the full-width-at-half-maximum (FWHM).

IXS directly probes the phonon dynamical structure factor,
$S(\mathbf{Q},\omega)$. The IXS cross-section for solid angle $d\Omega$ and
bandwidth $d\omega$ can be expressed as [41]:

$\frac{d^{2}\sigma}{d\Omega
d\omega}=\frac{k_{f}}{k_{i}}r_{0}^{2}|\vec{\epsilon_{i}}\cdot\vec{\epsilon_{f}}|^{2}S(\mathbf{Q},\omega)$
(19)

where k and $\epsilon$ represent the scattering vector and x-ray polarization
and $i$ and $f$ denote initial and final states. $r_{0}$ is the classical
radius of the electron. In a typical measurement, the energy transfer $\omega$
is much smaller than the incident photon energy (23.71 keV in our study).
Therefore, the term $\frac{k_{f}}{k_{i}}\sim 1$, and
$\frac{d^{2}\sigma}{d\Omega d\omega}\propto S(\mathbf{Q},\omega)$.

$S(\mathbf{Q},\omega)$ is related to the imaginary part of the dynamical
susceptibility, $\chi^{\prime\prime}\left(\mathbf{Q},\omega\right)$, through
the fluctuation-dissipation theorem [41]:

$S\left(\mathbf{Q},\omega\right)=\frac{1}{\pi}\frac{1}{(1-e^{-\omega/k_{B}T})}\chi\prime\prime(\mathbf{Q},\omega)$
(20)

Where $\chi\prime\prime(\mathbf{Q},\omega)$ can be described by the damped
harmonic oscillator form [41], which has antisymmetric Lorentzian lineshape:

$\chi^{\prime\prime}(\mathbf{Q},\omega)=\sum_{i}I_{i}[\frac{\Gamma_{i}}{(\omega-\omega_{\mathbf{Q},i})^{2}+\Gamma_{i}^{2}}-\frac{\Gamma_{i}}{(\omega+\omega_{\mathbf{Q},i})^{2}+\Gamma_{i}^{2}}]$
(21)

here $i$ indexes the different phonon peaks.

The phonon peak can be extracted by fitting the IXS spectrum at constant-
momentum transfer Q, using Eq. (20) and (21). Due to the finite experimental
resolution, the IXS intensity is a convolution of $S(\mathbf{Q},\omega)$ and
the instrumental resolution function, $R(\omega)$:

$I(\mathbf{Q},\omega)=S(\mathbf{Q},\omega)\otimes R(\omega)$ (22)

Here $R(\omega)$ was determined by fitting elastic peak near a Bragg condition
(Fig. S5).

The dashed curves in Fig. 2(e)-(g) of the main text represent fits of the
experimental data following Eq. (20) to (22). In all these fittings, the
intrinsic phonon linewidth is much smaller than the energy resolution $\Delta
E\sim$1.3 meV. Therefore, all peak-widths are essentially resolution-limited.
An additional resolution-limited elastic peak is included to account for the
elastic diffuse scattering. The extracted peak positions are shown in Fig.
3(a) of the main text.

Figure S5: IXS spectrum at Q = (6, -3, 0) r.l.u.. The orange dashed line is a
pseudo-voigt function fitting of the experimental data at the room
temperature. The extracted energy resolution using the full width half maximum
(FWHM) is 1.3meV. The error bars represent one standard deviation assuming
Poisson counting statistics.

##  Appendix C Temperature dependence of the in-plane phonon modes

Figure S6: Temperature dependence of the in-plane acoustic phonon modes.
(a)-(c) Raw IXS spectra of the transverse acoustic (TA) and longitudinal
acoustic (LA) modes at high ($\geq 250$ K) and low temperature ($\leq 30$ K).
The dashed lines are the fits to the data. While the TA mode displays a
softening effect towards low temperature, the LA mode presents no phonon
anomalies. (d) The relative peak shift of TA phonons in sample#1 and sample#2.

Figure S6 present the temperature dependent IXS results on the two in-plane
acoustic phonon modes i.e., the transverse acoustic mode (TA) and the
longitudinal acoustic mode (LA). Neither the TA mode nor the LA mode display a
drastic change over temperature, which demonstrates that the in-plane acoustic
phonons maintain the nearly divergent condition down below the Kitaev
temperature ($T_{K}\sim 100$ K). In the TA mode, phonon softening effects
towards low temperature can be clearly resolved in the raw IXS spectra taken
from two different samples [Fig. S6 (a)(b)], which was reported and discussed
in ref. [31], whereas the LA mode present no such phonon anomaly. This
observation corroborates with the result shown in the main text that the TA
mode is closer to the divergent condition compared to the LA mode.

##  Appendix D More on the out-of-plane transverse acoustic phonon branch

Figure S7: Bose-factor corrected IXS spectrum taken from Q=(6,-3,0) r.l.u.
($\Gamma_{1}$) to (6,-3,1) r.l.u. ($\Gamma_{2}$) at $T=300$ K.

Figure S7 shows the full dispersion of the out-of plane transverse acoustic
phonon (TAo). The TAo mode, which corresponds to a shear displacement between
adjacent $\alpha$-RuCl3 layers, present an extremely low energy band top at
only 3 meV. This is consistent with the stacking faults found in this layered
material [42, 43].

## References

  * Ronning  _et al._ [2006] F. Ronning, R. W. Hill, M. Sutherland, D. G. Hawthorn, M. A. Tanatar, J. Paglione, L. Taillefer, M. J. Graf, R. S. Perry, Y. Maeno, and A. P. Mackenzie, Thermal conductivity in the vicinity of the quantum critical end point in ${\mathrm{sr}}_{3}{\mathrm{ru}}_{2}{\mathrm{o}}_{7}$, Phys. Rev. Lett. 97, 067005 (2006).
  * Bourgeois-Hope  _et al._ [2016] P. Bourgeois-Hope, S. Chi, D. A. Bonn, R. Liang, W. N. Hardy, T. Wolf, C. Meingast, N. Doiron-Leyraud, and L. Taillefer, Thermal conductivity of the iron-based superconductor fese: Nodeless gap with a strong two-band character, Phys. Rev. Lett. 117, 097003 (2016).
  * Grissonnanche  _et al._ [2019] G. Grissonnanche, A. Legros, S. Badoux, E. Lefrançois, V. Zatko, M. Lizaire, F. Laliberté, A. Gourgout, J.-S. Zhou, S. Pyon, T. Takayama, H. Takagi, S. Ono, N. Doiron-Leyraud, and L. Taillefer, Giant thermal Hall conductivity in the pseudogap phase of cuprate superconductors, Nature 571, 376 (2019).
  * Samajdar  _et al._ [2019] R. Samajdar, M. S. Scheurer, S. Chatterjee, H. Guo, C. Xu, and S. Sachdev, Enhanced thermal Hall effect in the square-lattice Néel state, Nature Physics 15, 1290 (2019).
  * Han  _et al._ [2019] J. H. Han, J.-H. Park, and P. A. Lee, Consideration of thermal hall effect in undoped cuprates, Phys. Rev. B 99, 205157 (2019).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006), january Special Issue.
  * Read and Green [2000] N. Read and D. Green, Paired states of fermions in two dimensions with breaking of parity and time-reversal symmetries and the fractional quantum hall effect, Phys. Rev. B 61, 10267 (2000).
  * Nasu  _et al._ [2016] J. Nasu, J. Knolle, D. L. Kovrizhin, Y. Motome, and R. Moessner, Fermionic response from fractionalization in an insulating two-dimensional magnet, Nature Physics 12, 912 (2016).
  * Banerjee  _et al._ [2017a] M. Banerjee, M. Heiblum, A. Rosenblatt, Y. Oreg, D. E. Feldman, A. Stern, and V. Umansky, Observed quantization of anyonic heat flow, Nature 545, 75 (2017a).
  * Leahy  _et al._ [2017] I. A. Leahy, C. A. Pocs, P. E. Siegfried, D. Graf, S.-H. Do, K.-Y. Choi, B. Normand, and M. Lee, Anomalous thermal conductivity and magnetic torque response in the honeycomb magnet $\alpha\text{$-$}{\mathrm{rucl}}_{3}$, Phys. Rev. Lett. 118, 187203 (2017).
  * Banerjee  _et al._ [2018] M. Banerjee, M. Heiblum, V. Umansky, D. E. Feldman, Y. Oreg, and A. Stern, Observation of half-integer thermal Hall conductance, Nature 559, 205 (2018).
  * Hentrich  _et al._ [2018] R. Hentrich, A. U. B. Wolter, X. Zotos, W. Brenig, D. Nowak, A. Isaeva, T. Doert, A. Banerjee, P. Lampen-Kelley, D. G. Mandrus, S. E. Nagler, J. Sears, Y.-J. Kim, B. Büchner, and C. Hess, Unusual phonon heat transport in $\alpha\text{$-$}{\mathrm{rucl}}_{3}$: Strong spin-phonon scattering and field-induced spin gap, Phys. Rev. Lett. 120, 117204 (2018).
  * Kasahara  _et al._ [2018] Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Majorana quantization and half-integer thermal quantum hall effect in a kitaev spin liquid, Nature 559, 227 (2018).
  * Yamashita  _et al._ [2020] M. Yamashita, J. Gouchi, Y. Uwatoko, N. Kurita, and H. Tanaka, Sample dependence of half-integer quantized thermal hall effect in the kitaev spin-liquid candidate $\alpha\text{$-$}{\mathrm{rucl}}_{3}$, Phys. Rev. B 102, 220404 (2020).
  * Yokoi  _et al._ [2021] T. Yokoi, S. Ma, Y. Kasahara, S. Kasahara, T. Shibauchi, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, C. Hickey, S. Trebst, and Y. Matsuda, Half-integer quantized anomalous thermal Hall effect in the Kitaev material candidate $\alpha$-RuCl3, Science 373, 568 (2021).
  * Czajka  _et al._ [2021] P. Czajka, T. Gao, M. Hirschberger, P. Lampen-Kelley, A. Banerjee, J. Yan, D. G. Mandrus, S. E. Nagler, and N. P. Ong, Oscillations of the thermal conductivity in the spin-liquid state of $\alpha$-RuCl3, Nature Physics 10.1038/s41567-021-01243-x (2021).
  * Ye  _et al._ [2020] M. Ye, R. M. Fernandes, and N. B. Perkins, Phonon dynamics in the kitaev spin liquid, Phys. Rev. Research 2, 033180 (2020).
  * Feng  _et al._ [2021] K. Feng, M. Ye, and N. B. Perkins, Temperature evolution of the phonon dynamics in the kitaev spin liquid, Phys. Rev. B 103, 214416 (2021).
  * Takagi  _et al._ [2019] H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of kitaev quantum spin liquids, Nature Reviews Physics 1, 264 (2019).
  * Plumb  _et al._ [2014] K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, $\alpha-{\mathrm{rucl}}_{3}$: A spin-orbit assisted mott insulator on a honeycomb lattice, Phys. Rev. B 90, 041112 (2014).
  * Sandilands  _et al._ [2015] L. J. Sandilands, Y. Tian, K. W. Plumb, Y.-J. Kim, and K. S. Burch, Scattering continuum and possible fractionalized excitations in $\alpha\text{$-$}{\mathrm{rucl}}_{3}$, Phys. Rev. Lett. 114, 147201 (2015).
  * Banerjee  _et al._ [2016] A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Proximate kitaev quantum spin liquid behaviour in a honeycomb magnet, Nature Materials 15, 733 (2016).
  * Zhou  _et al._ [2016] X. Zhou, H. Li, J. A. Waugh, S. Parham, H.-S. Kim, J. A. Sears, A. Gomes, H.-Y. Kee, Y.-J. Kim, and D. S. Dessau, Angle-resolved photoemission study of the kitaev candidate $\alpha-{\mathrm{rucl}}_{3}$, Phys. Rev. B 94, 161106 (2016).
  * Banerjee  _et al._ [2017b] A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Neutron scattering in the proximate quantum spin liquid $\alpha$-rucl3, Science 356, 1055 (2017b).
  * Janša  _et al._ [2018] N. Janša, A. Zorko, M. Gomilšek, M. Pregelj, K. W. Krämer, D. Biner, A. Biffin, C. Rüegg, and M. Klanjšek, Observation of two types of fractional excitation in the kitaev honeycomb magnet, Nature Physics 14, 786 (2018).
  * Widmann  _et al._ [2019] S. Widmann, V. Tsurkan, D. A. Prishchenko, V. G. Mazurenko, A. A. Tsirlin, and A. Loidl, Thermodynamic evidence of fractionalized excitations in $\alpha-\mathrm{RuC}{\mathrm{l}}_{3}$, Phys. Rev. B 99, 094415 (2019).
  * Jackeli and Khaliullin [2009] G. Jackeli and G. Khaliullin, Mott insulators in the strong spin-orbit coupling limit: From heisenberg to a quantum compass and kitaev models, Phys. Rev. Lett. 102, 017205 (2009).
  * Laurell and Okamoto [2020] P. Laurell and S. Okamoto, Dynamical and thermal magnetic properties of the Kitaev spin liquid candidate $\alpha$-RuCl3, npj Quantum Materials 5, 2 (2020).
  * Sears  _et al._ [2020] J. A. Sears, L. E. Chern, S. Kim, P. J. Bereciartua, S. Francoual, Y. B. Kim, and Y.-J. Kim, Ferromagnetic Kitaev interaction and the origin of large magnetic anisotropy in $\alpha$-RuCl3, Nature Physics 16, 837 (2020).
  * Do  _et al._ [2017] S.-H. Do, S.-Y. Park, J. Yoshitake, J. Nasu, Y. Motome, Y. S. Kwon, D. T. Adroja, D. J. Voneshen, K. Kim, T.-H. Jang, J.-H. Park, K.-Y. Choi, and S. Ji, Majorana fermions in the Kitaev quantum spin system $\alpha$-RuCl3, Nature Physics 13, 1079 (2017).
  * Li  _et al._ [2021] H. Li, T. T. Zhang, A. Said, G. Fabbris, D. G. Mazzone, J. Q. Yan, D. Mandrus, G. B. Halász, S. Okamoto, S. Murakami, M. P. M. Dean, H. N. Lee, and H. Miao, Giant phonon anomalies in the proximate kitaev quantum spin liquid $\alpha$-rucl3, Nature Communications 12, 3513 (2021).
  * Bruin  _et al._ [2021] J. Bruin, R. Claus, Y. Matsumoto, N. Kurita, H. Tanaka, and H. Takagi, Robustness of the thermal hall effect close to half-quantization in a field-induced spin liquid state (2021), arXiv:2104.12184 [cond-mat.str-el] .
  * Lefrancois  _et al._ [2021] E. Lefrancois, G. Grissonnanche, J. Baglo, P. Lampen-Kelley, J. Yan, C. Balz, D. Mandrus, S. E. Nagler, S. Kim, Y.-J. Kim, N. Doiron-Leyraud, and L. Taillefer, Introduction to high-resolution inelastic x-ray scattering (2021), arXiv:2011.05493 [cond-mat.mtrl-sci] .
  * May  _et al._ [2020] A. F. May, J. Yan, and M. A. McGuire, A practical guide for crystal growth of van der waals layered materials, Journal of Applied Physics 128, 051101 (2020), https://doi.org/10.1063/5.0015971 .
  * [35] See Supplemental Materials at http…. .
  * Fåk and Dorner [1997] B. Fåk and B. Dorner, Phonon line shapes and excitation energies, Physica B: Condensed Matter 234, 1107 (1997).
  * Le Tacon  _et al._ [2011] M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. M. Sala, V. Hinkov, M. W. Haverkort, M. Minola, M. Bakr, K. J. Zhou, S. Blanco-Canosa, C. Monney, Y. T. Song, G. L. Sun, C. T. Lin, G. M. De Luca, M. Salluzzo, G. Khaliullin, T. Schmitt, L. Braicovich, and B. Keimer, Intense paramagnon excitations in a large family of high-temperature superconductors, Nat. Phys. 7, 725 (2011).
  * Miao  _et al._ [2018] H. Miao, D. Ishikawa, R. Heid, M. Le Tacon, G. Fabbris, D. Meyers, G. D. Gu, A. Q. R. Baron, and M. P. M. Dean, Incommensurate phonon anomaly and the nature of charge density waves in cuprates, Phys. Rev. X 8, 011008 (2018).
  * Weber  _et al._ [2011] F. Weber, S. Rosenkranz, J.-P. Castellan, R. Osborn, R. Hott, R. Heid, K.-P. Bohnen, T. Egami, A. H. Said, and D. Reznik, Extended phonon collapse and the origin of the charge-density wave in $2h\mathrm{\text{$-$}}{\mathrm{nbse}}_{2}$, Phys. Rev. Lett. 107, 107403 (2011).
  * Hoesch  _et al._ [2009] M. Hoesch, A. Bosak, D. Chernyshov, H. Berger, and M. Krisch, Giant kohn anomaly and the phase transition in charge density wave ${\mathrm{zrte}}_{3}$, Phys. Rev. Lett. 102, 086402 (2009).
  * Baron [2020] A. Q. R. Baron, Introduction to high-resolution inelastic x-ray scattering (2020), arXiv:1504.01098 [cond-mat.mtrl-sci] .
  * Sears [2017] J. Sears,  _Neutron and X-ray Diffraction Studies of Structure and Magnetism in $\mathrm{\alpha}{\mathrm{-}}{\mathrm{RuCl}}_{3}$_, Ph.D. thesis (2017).
  * Mi  _et al._ [2021] X. Mi, X. Wang, H. Gui, M. Pi, T. Zheng, K. Yang, Y. Gan, P. Wang, A. Li, A. Wang, L. Zhang, Y. Su, Y. Chai, and M. He, Stacking faults in $\mathrm{\alpha}{\mathrm{-}}{\mathrm{rucl}}_{3}$ revealed by local electric polarization, Phys. Rev. B 103, 174413 (2021).