# Robustness for Free: Quality-Diversity Driven Discovery of Agile Soft
Robotic Gaits

John Daly1, Daniel Casper1, Muhammad Farooq1, Andrew James1, Ali Khan1,
Phoenix Mulgrew1, Daniel Tyebkhan1, Bao Vo1, John Rieffel2 1Undergraduate
Student, Union College, Schenectady, NY 2 Professor, Computer Science
Department, Union College, Schenectady, NY<EMAIL_ADDRESS>
###### Abstract

Soft robotics aims to develop robots able to adapt their behavior across a
wide range of unstructured and unknown environments. A critical challenge of
soft robotic control is that nonlinear dynamics often result in complex
behaviors hard to model and predict. Typically behaviors for mobile soft
robots are discovered through empirical trial and error and hand-tuning. More
recently, optimization algorithms such as Genetic Algorithms (GA) have been
used to discover gaits, but these behaviors are often optimized for a single
environment or terrain, and can be brittle to unplanned changes to terrain. In
this paper we demonstrate how Quality Diversity Algorithms, which search of a
range of high-performing behaviors, can produce repertoires of gaits that are
robust to changing terrains. This robustness significantly out-performs that
of gaits produced by a single objective optimization algorithm.

###### Index Terms:

Soft Robotics, Quality Diversity Algorithms, Evolutionary Robotics

##  I Introduction

The field of soft robotics seeks to create bio-inspired soft systems whose
performance rivals that of natural systems. Such robots hold a particular
appeal for tasks in unstructured and uncertain environments, where compliance
and deformability might allow a robot to change shape to squeeze through small
apertures, or survive large drops. However, because they are nature high
dimensional dynamic systems with an essentially infinite number of degrees of
freedom, they are difficult to control by conventional means[1, 2].

As a result of this complexity, the gaits of mobile soft robots are typically
discovered through empirical trial-and-error[2], a process that can be both
challenging and time consuming. More recently, optimization algorithms such
have Genetic Algorithms have been used to automatically generate gaits for
soft robots [3, 4, 5, 6] but like the hand-tuned ones they tend to be
optimized for only a single terrain, and can therefore be brittle to changes.
This brittleness is a challenge to creating soft robots able to autonomously
adapt their behavior as needed, for instance when they encounter a novel
environmental substrate, or when they are physically damaged.

One promising alternative approach is to use novelty-oriented approaches such
as Quality Diversity Algorithms (QDAs) to develop repertoires of effective but
behaviorally distinct gaits. One advantage of these approaches might be that
the behavioral diversity of the repertoire offers some insurance to changes in
terrain. In this paper we demonstrate the non-robust nature of solutions
evolved with a state-of-the-art optimization algorithm (CMA-ES), and explore
how Quality Diversity Algorithms can not just discover effective gaits, but
can also be robust to significant changes in their environment.

##  II Background and Related Work

###  II-A Soft Robotics Simulators

Although the complexity of soft materials makes soft robots incredibly
difficult to simulate with fidelity sufficient for “sim-to-real transfer” [7]
one can still glean significant insights about soft robotics by generating
gaits in simulation. Several soft robotic simulators exist, including VoxCraft
[8, 7], SOFA [9, 10], EvoGym [6] and more recently differentiable Mass
Particle Method (MPM) simulators such as ChainQueen [11] and DiffTaiChi [12].

For this work we rely upon 2D-VSR-Sim, a 2D Java-based simulator that models a
single voxel as a deformable square with four rigid masses connected by both
spring–damper systems (SDSs) and tensile cords [13, 14]. Robot morphologies in
the simulator are constructed by connecting multiple voxels together at the
edges. The coefficients (mass, spring constants, etc) can be independently
configured for each voxel, allowing for a significant range of soft material
properties. Robots in 2D-VSR-Sim are typically controlled by applying
oscillating volumetric changes to individual voxels. These oscillations can
controlled in either an open-loop or closed-loop fashion during simulation.

###  II-B Optimization Methods

Because of the complexity of soft materials and the non-intuitive nature of
their behavior, optimization techniques such as Genetic Algorithms (GAs) are a
promising approach for the automated discovery of soft robotic gaits. One
advantage is that they can reduce the bias inherent in the human-oriented
empirical trial-and-error that is often used in soft robotics. GAs have been
successfully used to evolve soft robotic gaits [3, 4, 5, 6]. Because they
often require thousands (or more) evaluations to discover effective behaviors,
GAs and other optimizers often rely simulators to produce results.

For this work we consider two types of optimization method:

#### CMA-ES

: Covariance Matrix Adaptation Evolution Strategies (CMA-ES) is a well
established off-the-shelf single-objective optimizer. We use it in our work as
a baseline optimizer against which to compare our QDA approach.

#### Quality Diversity Algorithms

Unlike single-objective optimizers like CMA-ES, Quality Diversity algorithms
(QDAs) [15, 16, 17, 18] search for a combination of both novelty and quality.
When applied to robotics, this means finding mappings between a robot’s
parameter space (those features that can be affected by a controller), and its
behavior space (a description of the outcomes of those actions). The outcome
of a quality diversity search is an archive describing the variety of
behaviors that the robot is capable of – essentially a mapping between
parameter space and behavior space.

QD algorithms have been employed with considerable success with both simulated
virtual agents [17, 16] and physical robots [18]. Most notably, Cully et al.
[18] used a hybrid simulation/physical approach to behavior build repertoires
for a hexapod robot.

Recently, in our own work [19] we explored methods to autonomously discover
novel and effective soft robotic locomotive behaviors by using QDAs.
Specifically, we demonstrated the ability to to autonomously discover a
diverse repertoire of unique locomotive behaviors for a tether-free soft
tensegrity robot. The behavior space in question described the translational
and rotational displacement of the robot along the 2D plane. These discovered
behaviors could, in turn, be incorporated into higher-level control strategies
such as policy networks [16].

Our interest in this work is not about discovering a repertoire of behaviors
producing different displacements, but instead to explore the variety of
qualitatively different and interesting gaits a 2D voxel-based robot can
traverse a horizontal landscape.

##  III Methods

Our ambition was to explore the robustness of gaits optimized via a
conventional single-objective optimizer (CMA-ES) against the archives of
behaviors produced by a Quality Diversity Algorithm, specifically the popular
QDpy package[20]. Robustness in this case means the loss (or gain) of fitness
when a gait evolved on one terrain is transferred to a novel terrain without
any subequent re-training.

To explore this, we ran each optimization method on a fixed robot morphology
across a variety of terrains, and explored how each optimized system responded
to subsequent post-optimization changes in terrain. In each case, the
optimization method was given a fixed evaluation budget in order to make
meaningful comparisons between the two approaches.

###  III-A Designing Terrains

The first step for designing this experiment was creating unique terrains for
the gaits to be evolved on. Terrains were generated as periodic waves, and
described in 2D-VSR-Sim as a sequence of coordinate points along the plane.
Terrains were chosen to be unique and variable, but not impossible to
traverse. The surfaces were thefore smaller in dimension than the robot
itself. The six terrains chosen, shown in Figure 1 were:

  * •

flat: a flat horizontal surface

  * •

spiky: periodic triangles 1 unit wide and 0.5 units tall

  * •

longspikes: periodic triangles 2 units wide and 0.5 units tall.

  * •

longerspikes: periodic triangles 4 units wide and 0.5 units tall.

  * •

sparsespike: a periodic asymmetric sawtooth shape 1.5 units wide and 0.5 units
tall.

  * •

valley: an aperiodic valley sloping upwards from the origin in each direction.

(a) Flat

(b) Spiky

(c) Longspikes

(d) Longerspikes

(e) Sawtooth

(f) Valley

Figure 1: The six terrains used in this study.

###  III-B Robot Construction and Actuation

As described above, 2D-VSR-Sim is a voxel-based simulator, and robot
morphologies consist of groupings of connected voxels which can be
individually actuated (or passive). For a given robot design, actuation is
achieved by expanding and contracting the actuated voxels. Details on the
simulator are provided in [13].

The robots used in this experiment were constructed from 5 voxels in a simple
biped-shape, as shown in Figure 2

Figure 2: The anatomy and shape of the 5-voxel robot used in this work. Each
of the three column of voxels for the robot was independently actuated with a
sinusoidal wave.

For this work the robot was actuated by a sinusoidal oscillators controlling
each of the three vertical columns of voxels. The specific actuation
parameters varied in order to produce gaits were:

  * •

amplitude: a global property shared by all voxels describing how much voxels
expand and contract over a given oscillation.

  * •

frequency: a second global value describing the rate at which all voxels
oscillate.

  * •

column phase: each vertical column of the robot has its own phase offset,
causing different columns of the robot to oscillate differently.

These five parameters, encoded as floating point values in the range (0,1),
were used both as the genotype for CMA-ES and the parameter space for QDPy.

###  III-C Fitness

The fitness, used for both CMA-ES and QDA was the speed (in units per second)
at which the robot traveled horizontally during a gait’s simulated evaluation
period. Absolute horizontal speed was used here to avoid biasing solutions
towards a “front” or “back” of the robot – having a behavior that moves in a
subjectively “backwards” direction, for instance to the left, is simply a
matter of perspective.

###  III-D Behavioral Descriptors

The one feature that distinguishes QDAs from optimizers is the use of a
behavior-space descriptor, which maps some qualitative measure of how a gait
behaves down into a vector of numerical values. These descriptor vectors are
then used to populate the archive of elite solutions.

In order to facilitate searching for quality alongside diversity though, it’s
necessary to represent the behavior space as a discrete mapping of “behavioral
niches” or bins in which all behaviors inside the niche are similar enough
that one behavior can reasonably represent the whole subset. The behavior that
lies within a niche and is the best performing seen so far is called the elite
of the niche. The collection of all the elites is known as an archive. Our
behavior space for this work is divided into 10 regions along each axis,
creating 100 distinct niches, as illustrated by Figure 4.

Because our interest was in describing different properties of a soft robotic
gait we chose descriptors that measured various characteristics of soft
bodies. These descriptors, as illustrated in Figure 3 are as follows:

  * •

squish: which measures the variance in the distance between the bottom-left
(south-west) and top-right (north-east) corners of the robot. A gait that
exhibits low “squish” has low variance in this distance during evaluation,
whereas a robot with high “squish” varies significantly.

  * •

wobble: which measures the variance in rotation, or pitch of the robot
relative to the center. A low-wobble gait remains relatively steady in pitch
during evaluation, whereas a high-wobble gait tips back and forth
significantly.

Figure 3: Two qualitative measures were used to describe a gait, “squish”
(variance in distance between two corners), and “wobble” (variance in rotation
around the center of mass. These descriptors were used to place gaits into the
QDA archive. Figure 4: The Quality Diversity Algorithm generates an archive of
high fitness but behaviorally distinct. The motion of the robot during the
gait is classified along two dimensions - “squish” (variance in diagonal
width) and “wobble” (variation in pitch around the center of mass of the
robot), and then placed into a corresponding discreted bin. Whenever two gaits
are placed into the same bin of the archive only the fittest of the two is
retained in the archive.

###  III-E Experimental Procedure

With the parameters determining actuation established, as well as a
quantitative measure of fitness and qualitative measures of behavior, we can
now compare the robustness of gaits evolved under the two treatments: CMA-ES
and QDA. Specifically we are seeking to measure how gaits evolved on one
terrain perform when transferred to a novel terrain.

Each algorithm was given the same budget for evolution: 25 seconds of
simulation time for each individual evaluation, and a maximum of 2000
evaluations during each optimizing run. Each algorithm was also given the same
population/batch size (20).

In the case of CMA-ES this means that the optimizer ran with a population size
of 20 and for a maximum of 2000 total evaluations, using horizontal speed as a
fitness function.

In the case of the QDA, this means that algorithm ran with a batch size of 20
for a maximum of 2000 total evaluations, using horizontal speed as fitness,
and the two behavioral descriptors described above to define the archive.

We then followed these steps to compare the robustness across terrains of the
QDA and CMA-ES, as described in the algorithm below. The process can be broken
down into th phases. For an optimization phase, each algorithm was run for its
complete evaluation budget on a pre-determined terrain, recording the highest
fitness gait (CMA-ES) and highest-fitness gait from the archive (QDA) at the
end. For a transfer phase, the highest-fitness CMA-ES gait was re-evaluated on
each of the other terrains, recording the new fitness. By contrast, the entire
archive of elites from the QDA was re-evaluated, and the fitness of the
highest-performing elite was recorded. We repeated this 30 times for each
terrain, and repeated the entire process across all six terrains.

procedure Experimental Procedure

$T\leftarrow All\>Defined\>Terrains$

for $t\in T$ do

(Optimization Phase)

for $n=1,\dots,30trials$ do

Generate QDA gait archive on terrain $t$

Evolve CMA-ES gait on terrain $t$

(Transfer Phase)

$R\leftarrow T-t$

for $r\in R$ do

Retest collected QDA archives on $r$

Retest collected CMA-ES behaviors on $r$

All experiments were run on baremetal compute instances provided by the
Chameleon Cloud project [21].

With these data gathered we can now compare the fitness gain/loss for each
CMA-ES gait before and after transfer against the best new fitness from the
corresponding QDA elite archive, averaged across all 30 trials.

##  IV Results

Figures 5 and 6 illustrate the results of this experiment. Let us begin with
Figure 5. The left hand column of each subfigure shows the maximum fitness
achieved by each algorithm on the initial terrain. Not surprisingly, in every
case CMA-ES substantially outperforms the Quality Diversity Algorithm. This is
intuitive, because the CMA-ES is able to focus its entire evaluation budget
myopically on optimizing fitness across that one terrain, whereas the QDA has
to spend its evaluation budget creating an entire archive of fit but
behaviorally distinct solutions.

The remaining columns in each subfigure demonstrate the loss/gain of fitness
when the solutions from the optimization phase are transferred to each of the
novel terrains. Across almost all treatments, the CMA-ES fitness drops
substantially in the novel terrain, illustrating the brittle and non-robust
nature of these evolved solutions. The exception in this case is the “valley”
terrain, which is substantially more difficult than all the other terrains -
meaning that the CMA-ES solutiosn improve on the other terrains. By contrast,
the best fitness from the QDA elite archive varies substantially less - and in
many cases increases on the novel terrain. There are notably increases in QDA
elite fitness on new terrains - most often when QDA elites evolved on a non-
flat terrain are transferred to the “flat” terrain.

Figure 6 shows another perspective on these same data, aggregating the fitness
gain/loss for each initial terrain. For CMA-ES (left) we can see that, with
the exception of those evolved on the “valley” terrain all evolved solutions
are not robust when transferred to a novel terrain. Again, however for the QDA
archives (right), the fitness losses are less extreme. In fact, unlike for the
CMA-ES solutions, all of the QDA archives perform better when transferred to
the flat terrain. More interestingly the “spiky” terrain seems to generate
archives that are more robust to novel terrains than the other terrains.

Figure 5: Results of generating solutions one one terrain (left-most column in
each subfigure and then re-evaluating those solutions on each of the other
terrains. While CMA-ES gaits generally outperforms the QDA, they are quite
brittle when transferred to new terrains. QDA results have lower initial
fitness, but generally are much more robust to changes in terrain.

Figure 6: Comparisons of absolute gain (blue) or loss (red) of fitness when
gaits generated on one terrain are transferred to each other terrain. The QDA
solutions (right) lose less - and gain more - fitness when transferred than
the CMA-ES solutions

##  V Discussion and Conclusion

In this work we’ve compared the robustness soft robotic gaits evolved with a
single-objective optimizer, CMA-ES, against the archives of gaits generated by
a Quality Diversity Algorithm. When given the same budget of gait evaluations
CMA-ES significantly outperforms QDAs in each of the terrains.

However when those CMA-ES solutions are transferred to novel terrains their
fitness drops significantly. This in of itself is not a surprising result:
CMA-ES solutions are specialized to work on the terrain they were trained on.
The archives of diverse solutions generated by the Quality Diversity
Algorithm, by contrast, are more robust to transfer. This does not mean that
the archive elite that was best on the original terrain also does well on the
novel terrain - instead it means one qualitatively distinct gaits in the
archive does comparatively well on the novel terrain. In other words, the
diversity of “wobble” and “squish” in the archive allows for enough some
benefit when transferred to a new terrain.

The one exception to this pattern is the “valley” terrain - one that both CMA-
ES and QDA approaches struggle to find effective gaits for. Nonetheless, the
CMA-ES solutions actually increase in fitness when transferred to other
terrains. The QDA solutions generated on the spiky terrain by contrast only
improves on the flat terrain.

Compellingly, the “spiky” terrain seems to elicit the most robustness from the
QDA archives, which drop in fitness only on the challenging “valley” terrain.

We can draw two insights from these results. The first is that the time spent
generating diverse archives of gaits with QDAs has on one hand a cost in terms
of discovering less optimal solutions than CMA-ES. However this investment in
diversity pays dividends in the long run, by providing some inherent
robustness to terrain change. This could be leveraged in future work, for
instance, by using a QDA to develop an archive of behaviorally distinct gaits
for a soft robot, and using a policy network to determine which gait is most
appropriate when a terrain changes. The second insight is that the choice of
terrain that gaits are generated in can have a significant effect upon the
ability of those solutions to robustly transfer to a new terrain. For CMA-ES
this happens to the the challenging “valley” terrain, and for QDA it is the
spiky terrain. In future work we hope to explore how to find terrains that can
somehow maximize the robustness of a generated archive.

## Acknowledgment

This work was supported by NSF Award ID NRI-1939930, and by the Union College
Summer Research Fellowhsip program. The authors would also like to thank the
ChameleonCloud project for access to reliable high-speed baremetal compute
instances.

## References

  * [1] H. Lipson, “Challenges and opportunities for design, simulation, and fabrication of soft robots,” _Soft Robotics_ , vol. 1, no. 1, pp. 21–27, 2014.
  * [2] R. F. Shepherd, F. Ilievski, W. Choi, S. A. Morin, A. A. Stokes, A. D. Mazzeo, X. Chen, M. Wang, and G. M. Whitesides, “Multigait soft robot,” _Proceedings of the National Academy of Sciences_ , vol. 108, no. 51, pp. 20 400–20 403, 2011.
  * [3] J. Rieffel, F. Saunders, S. Nadimpalli, H. Zhou, S. Hassoun, J. Rife, and B. Trimmer, “Evolving soft robotic locomotion in physx,” in _Proceedings of the 11th annual conference companion on genetic and evolutionary computation conference: Late breaking papers_ , 2009, pp. 2499–2504.
  * [4] J. Rieffel, D. Knox, S. Smith, and B. Trimmer, “Growing and evolving soft robots,” _Artificial life_ , vol. 20, no. 1, pp. 143–162, 2014.
  * [5] F. Cochevelou, D. Bonner, and M.-P. Schmidt, “Differentiable soft-robot generation,” in _Proceedings of the Genetic and Evolutionary Computation Conference_ , 2023, pp. 129–137.
  * [6] J. Bhatia, H. Jackson, Y. Tian, J. Xu, and W. Matusik, “Evolution gym: A large-scale benchmark for evolving soft robots,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 2201–2214, 2021.
  * [7] S. Kriegman, A. M. Nasab, D. Shah, H. Steele, G. Branin, M. Levin, J. Bongard, and R. Kramer-Bottiglio, “Scalable sim-to-real transfer of soft robot designs,” in _2020 3rd IEEE International Conference on Soft Robotics (RoboSoft)_. IEEE, 2020, pp. 359–366.
  * [8] S. Liu, D. Matthews, S. Kriegman, and J. Bongard, “Voxcraft-sim, a gpu-accelerated voxel-based physics engine,” https://github.com/voxcraft/voxcraft-sim, 2020.
  * [9] C. Duriez, “Control of elastic soft robots based on real-time finite element method,” in _2013 IEEE international conference on robotics and automation_. IEEE, 2013, pp. 3982–3987.
  * [10] S. E. Navarro, S. Nagels, H. Alagi, L.-M. Faller, O. Goury, T. Morales-Bieze, H. Zangl, B. Hein, R. Ramakers, W. Deferme _et al._ , “A model-based sensor fusion approach for force and shape estimation in soft robotics,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5621–5628, 2020.
  * [11] Y. Hu, J. Liu, A. Spielberg, J. B. Tenenbaum, W. T. Freeman, J. Wu, D. Rus, and W. Matusik, “Chainqueen: A real-time differentiable physical simulator for soft robotics,” in _2019 International conference on robotics and automation (ICRA)_. IEEE, 2019, pp. 6265–6271.
  * [12] Y. Hu, L. Anderson, T.-M. Li, Q. Sun, N. Carr, J. Ragan-Kelley, and F. Durand, “Difftaichi: Differentiable programming for physical simulation,” _arXiv preprint arXiv:1910.00935_ , 2019.
  * [13] E. Medvet, A. Bartoli, A. De Lorenzo, and S. Seriani, “2d-vsr-sim: A simulation tool for the optimization of 2-d voxel-based soft robots,” _SoftwareX_ , vol. 12, p. 100573, 2020.
  * [14] ——, “Design, validation, and case studies of 2d-vsr-sim, an optimization-friendly simulator of 2-d voxel-based soft robots,” _arXiv preprint arXiv:2001.08617_ , 2020.
  * [15] J. K. Pugh, L. B. Soros, and K. O. Stanley, “Quality diversity: A new frontier for evolutionary computation,” _Frontiers in Robotics and AI_ , vol. 3, p. 40, 2016. [Online]. Available: https://www.frontiersin.org/article/10.3389/frobt.2016.00040
  * [16] M. Duarte, J. Gomes, S. M. Oliveira, and A. L. Christensen, “Evolution of repertoire-based control for robots with complex locomotor systems,” _IEEE Transactions on Evolutionary Computation_ , vol. 22, no. 2, pp. 314–328, 2018.
  * [17] J. Lehman and K. O. Stanley, “Evolving a diversity of virtual creatures through novelty search and local competition,” in _Proceedings of the 13th annual conference on Genetic and evolutionary computation_. ACM, 2011, pp. 211–218.
  * [18] A. Cully, J. Clune, D. Tarapore, and J.-B. Mouret, “Robots that can adapt like animals,” _Nature_ , vol. 521, no. 7553, pp. 503–507, 2015.
  * [19] K. Doney, A. Petridou, J. Karaul, A. Khan, G. Liu, and J. Rieffel, “Behavioral repertoires for soft tensegrity robots,” in _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_. IEEE, 2020, pp. 2265–2271.
  * [20] L. Cazenille, “Qdpy: A python framework for quality-diversity,” https://gitlab.com/leo.cazenille/qdpy, 2018.
  * [21] K. Keahey, J. Anderson, Z. Zhen, P. Riteau, P. Ruth, D. Stanzione, M. Cevik, J. Colleran, H. S. Gunawi, C. Hammock _et al._ , “Lessons learned from the chameleon testbed,” in _2020 USENIX annual technical conference (USENIX ATC 20)_ , 2020, pp. 219–233.


Signatures of fractionalization in the optical phonons of hyperhoneycomb  Kitaev  magnet β-Li_2IrO_3
Swetlana Swarup
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Susmita Singh
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Natalia B. Perkins
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

In this study, we propose that the signatures of spin fractionalization in quantum magnets can be identified through a detailed analysis of the temperature dependence of the asymmetric Fano lineshape of optical phonons overlapping with a continuum of spin excitations. We focus on the hyperhoneycomb magnet $\beta$-Li$_2$IrO$_3$, a promising candidate for being in proximity to a three-dimensional Kitaev quantum spin liquid. The Raman response in $\beta$-Li$_2$IrO$_3$ notably displays a distinctive asymmetric Fano lineshape in the 24 meV Raman-active optical phonon. This asymmetry arises from the interaction between the discrete phonon mode and the spin excitation continuum, which could be fractionalized if the material is indeed near a quantum spin-liquid phase. Our theoretical model considers the coupling of this optical phonon to Majorana fermions in the Kitaev model on the hyperhoneycomb lattice. Our findings reveal that the temperature-dependent Fano lineshape is consistent with the fractionalization of spins into Majorana fermions and $\mathbb{Z}{_2}$ fluxes.

§ INTRODUCTION

Quantum spin liquids (QSLs) are among the most fascinating phases of matter due to their fractionalized quasiparticles and nontrivial topological properties <cit.>.
Among the various types of QSLs, the exactly solvable Kitaev model has garnered significant attention for its potential realization in real materials
<cit.>. Numerous candidate materials have been proposed, including
$\alpha-\text{RuCl}_3$, $\text{Na}_2\text{IrO}_3, \text{Cu}_3\text{NaIr}_2\text{O}_6$ and $\alpha$-, $\beta$-, $\gamma$-$\text{Li}_2\text{IrO}_3$, $\text{H}_3\text{LiIr}_2\text{O}_6,$
to name a few <cit.>.

However, identifying  any  QSLs, including Kitaev QSLs, in real materials poses a unique challenge because they are difficult to detect through direct experimental probes <cit.>. Nevertheless, progress has been made by combining information from various probes such as specific heat measurements, inelastic neutron and Raman scattering, resonant inelastic X-ray scattering (RIXS), and thermal transport to gain insights into the different properties that characterize a QSL <cit.>.

Recently, the study of phonon dynamics has emerged as a powerful tool to uncover the dynamics of underlying fractionalized degrees of freedom <cit.>. In particular, it was observed that sound attenuation measured in $\alpha$-RuCl$_3$ through ultrasound experiments <cit.> shows characteristic temperature behavior that supports the scattering of acoustic phonons from Majorana fermions in the extended Kitaev spin liquid, as predicted in earlier theoretical studies <cit.>. Additionally, both theoretical <cit.>  and experimental studies <cit.> have noted that optical phonons in QSL candidate materials, probed through Raman spectroscopy, exhibit a characteristic Fano lineshape. This lineshape arises from the overlap of phonons with an underlying continuum of fractionalized excitations. Its evolution with temperature and in the presence of an external magnetic field can provide valuable insights into the nature of these fractionalized excitations.

Motivated by these studies, in this paper
we investigate the low-energy optical phonon dynamics
in  $\beta$-Li$_2$IrO$_3$ <cit.>, a promising platform for the relatively rare QSL behavior on a three-dimensional lattice with bond-anisotropic Ising-like  interactions <cit.>.
The Kitaev model on the hyperhoneycomb lattice, describing these interactions,  is exactly solvable, featuring a QSL ground state with fractionalized excitations: gapless Majorana fermions and gapped $\mathbb{Z}_2$ fluxes. The ground state is in the zero-flux sector, with a finite-temperature transition separating it from a high-temperature disordered flux state <cit.>. Notably, the Majorana fermions exhibit a nodal line structure <cit.>.
These inherent differences  from a two-dimensional counterpart
lead to fundamentally different thermodynamic behavior, especially at low temperatures.

Experimentally, a Fano lineshape has indeed been  observed in one of the low-energy phonons in $\beta$-Li$_2$IrO$_3$ <cit.>.
Despite the fact that the minimal spin Hamiltonian for $\beta$-Li$_2$IrO$_3$ includes non-Kitaev interactions in addition to the Kitaev coupling  <cit.>, leading to its
complex incommensurate
order with counter-rotating spirals
<cit.>, the material exhibits notable characteristics
of  a `proximate spin-liquid' regime with signatures of long-lived fractionalized excitations, which separate the low-temperature ordered phase from the high-temperature paramagnetic regime <cit.>.
For example, recent RIXS measurements <cit.> have shown that $\beta$-Li$_2$IrO$_3$ exhibits a characteristic
continuum response extending up to at least 300 K, well above its ordering temperature of 38 K. This supports theoretical predictions <cit.> and allows for the formation of the Fano lineshape in low-energy phonons <cit.>.

In this work, we theoretically investigate the temperature dependence of the parameters characterizing the Fano lineshape of a  24 meV phonon observed in Raman spectroscopy  in  $\beta$-Li$_2$IrO$_3$ <cit.>. To this end, we consider a spin-phonon Kitaev model on the hyperhoneycomb lattice and  systematically compute the Raman intensity.  Our focus is on Majorana fermionic degrees of freedom while we disregard the scattering processes from $\mathbb{Z}_2$ fluxes.  We show  that the spin-phonon coupling renormalizes phonon propagators and generates the salient Fano linshape observed in $\beta$-Li$_2$IrO$_3$.
Our numerical results for its temperature evolution give a good account of the  experimental observations  <cit.>, suggesting the proximity of $\beta$-Li$_2$IrO$_3$ to a QSL phase.

The contents of the paper are organized as follows: We introduce the spin-lattice coupled Kitaev Model in Sec. <ref> with the details of the spin, phonon and spin-phonon coupled Hamiltonian presented in Sec. <ref>, <ref> and <ref> respectively.  One-loop corrections to the phonon self-energy
due to spin-phonon interaction are calculated in Sec. <ref>.
In Sec. <ref>, we employ diagrammatic perturbation theory to evaluate the Raman intensity. We first obtain the pure magnetic Raman operator in Sec. <ref> and phonon Raman operator in Sec. <ref>. Then, in  Sec. <ref>,  we systematically calculate the Raman intensity in the spin-phonon coupled system using the S-matrix expansion.
We discuss the results in Sec. <ref>, where
we first determine the parameters of our theoretical model by fitting it to experimental data in Sec. <ref>. Next, we compute the temperature evolution of the Fano parameters in Sec. <ref>. Finally, we examine the dependence of these parameters on the Kitaev interaction strength in
Sec. <ref>. We conclude by summarizing our findings in Sec. <ref>.

§ MODEL

We consider the Kitaev model on the three-dimensional hyperhoneycomb lattice, incorporating coupling to low-energy optical phonons via the magnetoelastic interaction:

\begin{equation}
\mathcal{H}=\mathcal{H}_{\text{s}}+\mathcal{H}_{\text{ph}}+\mathcal{H}_{\text{s-ph}}
\end{equation}

Here, $\mathcal{H}_{\text{s}}$ represents the spin Hamiltonian, $\mathcal{H}_{\text{ph}}$ describes the phonon Hamiltonian, and $\mathcal{H}_{\text{s-ph}}$ captures the spin-phonon coupling.

§.§ Spin Hamiltonian

The spin Hamiltonian is given by the isotropic Kitaev model on the hyperhoneycomb lattice:

\begin{equation}\label{Hspin}
\mathcal{H}_{\text{s}}=-J_K \sum_t\sum_{\braket{ij}\in t} \sigma^{\alpha_t}_i\sigma^{\alpha_t}_j,\quad t \in \{x,x',y,y',z\}
\end{equation}

where $t$ represents the various nearest neighbor bonds connecting sites $i\in$ odd-sublattices and $j\in$ even-sublattices as shown in Fig. <ref> and $\alpha_t$ is the Cartesian component $\alpha_x=\alpha_{x'}=x$,   $\alpha_y=\alpha_{y'}=y$ and $\alpha_z=z$ associated with $t$ <cit.>.

Using Kitaev's representation of spins  expressed by the Pauli matrices $\boldsymbol{\sigma}=(\sigma^x,\sigma^y,\sigma^z)$
in terms of Majorana fermions, $\sigma_j^{\alpha_t}=ib_j^{\alpha_t} c_j$, and defining link variables as $\hat{u}_{ij}=ib_i^{\alpha_t}b_j^{\alpha_t}$, we obtain the free Majorana fermion Hamiltonian as

\begin{equation}\label{H-spin-free}
\mathcal{H}_{\text{s}}=i J_K\sum_t\sum_{\braket{ij}\in t} \, \hat{u}_{ij}\,c_i c_j.
\end{equation}

We restrict to the
zero-flux sector by setting all the link variables to +1.
Then, the free Hamiltonian (<ref>) can be easily diagonalized in Fourier space. Using the translational symmetry of the Bravais lattice and denoting the lattice point at site $i$ as ${\bf r}$, the Fourier transform is defined over the primitive unit cell as $c_{\mathbf{r},\delta}=\sqrt{\frac{2}{N}}\sum_\mathbf{k}\,e^{i\mathbf{k}\cdot\mathbf{r}}c_{\mathbf{k},\delta}$ and inversely, $c_{\mathbf{k},\delta}=\sqrt{\frac{1}{2N}}\sum_\mathbf{r}\,e^{-i\mathbf{k}\cdot\mathbf{r}}c_{\mathbf{r},\delta}$, where $\mathbf{r}=n_1\mathbf{a_1}+n_2\mathbf{a_2}+n_3\mathbf{a_3}$ and $\mathbf{k}=q_1\mathbf{k_1}+q_2\mathbf{k_2}+q_3\mathbf{k_3}$ for $n,q\in \mathbb{Z}$ and $N$ is the number of primitive unit cells. Here $\delta=1,2,3,4$ is the sublattice index in a given primitive unit cell.

The free Majorana fermion Hamiltonian in momentum-space takes the form

\begin{equation}
\mathcal{H}_{\text{s}}=\sum_\mathbf{k} C^\dagger_\mathbf{k}H_\mathbf{k}^{s}C_\mathbf{k}
\end{equation}

where $C_{\bf k}^\dagger=C_{\bf -k}$  with $C_{\bf k}^\dagger=(c_{\mathbf{-k},4},\, c_{\mathbf{-k},2}, \,c_{\mathbf{-k},1}, \, c_{\mathbf{-k},3})$. This arrangement is chosen such that the  matrix $H_\mathbf{k}^{s}$ takes a block diagonal form (<ref>).
To diagonalize this Hamiltonian, we express the Majorana fermions in terms of the complex fermions as

\begin{equation}\label{eqn:c2a}
C_{\bf k}=\sqrt{2} U_{\bf k} A_{\bf k}
\end{equation}

where  $U_\mathbf{k}$ is a unitary matrix that diagonalizes the matrix $H_{\bf k}^s$ and $A^\dagger_\mathbf{k}=(b_{-\mathbf{k}},\, a_{-\mathbf{k}},\, a_\mathbf{k}^\dagger ,\, b_\mathbf{k}^\dagger)$ where $a,b$ are the two flavors of the Bogoliubov complex fermions.

The spin Hamiltonian is then diagonal in the Bogoliubov basis and can be expressed as

\begin{equation}\label{eqn:spinhambogoliu}
\mathcal{H}_{\text{s}}=\sum_\mathbf{k} A^\dagger_\mathbf{k} \Tilde{H}^{\text{s}}_\mathbf{k} A_\mathbf{k},
\end{equation}

such that $\Tilde{H}^{\text{s}}_\mathbf{k}=2U_\mathbf{k}^\dagger H^{\text{s}}_{\mathbf{k}} U_\mathbf{k}$.  The eigenvalues are denoted by $\pm\epsilon_{\bf k}^a,\pm\epsilon_{\bf k}^b$, where $b$ represents the higher band and $a$ is the lower band as shown in Fig. <ref>. Further details on the geometry of the hyperhoneycomb lattice and the explicit form of the spin Hamiltonian can be found in Appendix <ref>.

Nearest neighbor bonds in the hyperhoneycomb lattice. A primitive unit cell (denoted in opaque colors) consists of 4 sublattices labeled $1-4$. Four primitive unit cells make up the orthorhombic unit cell. In orthorhombic coordinates (shown in red), the bonds are $\mathbf{M_z}=(0,0,1)$, $\mathbf{M_y}=(-\frac{1}{2},\frac{1}{\sqrt{2}},-\frac{1}{2})$, $\mathbf{M_{x'}}=(\frac{1}{2},-\frac{1}{\sqrt{2}},-\frac{1}{2})$, $\mathbf{M_{y'}}=(\frac{1}{2},\frac{1}{\sqrt{2}},-\frac{1}{2})$ and $\mathbf{M_x}=(-\frac{1}{2},-\frac{1}{\sqrt{2}},-\frac{1}{2})$. The crystallographic orthorhombic coordinates are related to the Cartesian axes (shown in blue) by $\mathbf{\hat{x}}=(\mathbf{\hat{a}+\hat{c}})/\sqrt{2},\mathbf{\hat{y}}=(\mathbf{\hat{c}-\hat{a}})/\sqrt{2}$ and $\mathbf{\hat{z}=-\hat{b}}$

§.§ Optical phonons

The free phonon Hamiltonian $\mathcal{H}_{\mathrm{ph}}\left(\mathcal{p}_i(\mathbf{r}),\mathcal{q}_i(\mathbf{r})\right)$ is a function of the canonical coordinate $\mathcal{q}_i(\mathbf{r})$,
which represents the lattice degrees of freedom in a unit cell at position ${\bf r}$,  and the corresponding canonical momentum $\mathcal{p}_i(\mathbf{r})$. In $\beta\text{-}\text{Li}_2\text{IrO}_3$, the unit cell contains
4 Iridium, 8 Lithium, and 12 Oxygen atoms.   Therefore, the canonical coordinate can be expressed as
$\mathbcal{q}(\mathbf{r})= (x_1, y_1, z_1, \ldots, x_{24}, y_{24}, z_{24})_{\bf r}$, representing the positions of the 24 atoms in the unit cell. We omit the explicit $\mathbf{r}$ dependence in the phonon displacement field because the long wavelength of the incident light leads to uniform lattice vibrations.

We are particularly interested in the optical phonons observed in the Raman spectra at zero momentum transfer, i.e., $\mathbf{q}=0$, originating from the center of the Brillouin zone (as shown
in Fig. <ref>). These modes can be described as a linear superposition of the displacement fields $u_{\Gamma}$, given by $u_{\Gamma } = \sum_{i=1}^{72} u_{\Gamma, i} \mathcal{q}_i$.  These displacement fields transform according to the irreducible representations (irreps) $\Gamma$ of the point group $D_{2h}$.
Focusing  on the 3D vibrations of the Ir atoms only, we restrict the sum to $i=12$. The  diagonalized Hamiltonian for $\mathbf{q}=0$ phonons  with energies $\Omega_\Gamma$ can now be written as:

\begin{equation}
\mathcal{H}_{\text{ph}}=\sum_{\Gamma}\Omega_\Gamma \bigg(\beta^\dagger_{\Gamma}\beta_{\Gamma}+\frac{1}{2}\bigg),
\end{equation}

where $\Gamma$ labels the irreps and $\beta^\dagger_{\Gamma}(\beta_{\Gamma})$ are the phonon creation (annihilation) operators, which are related to the displacement fields as

\begin{equation}
u_\Gamma(t)=i\bigg(\frac{\hbar}{2\rho \delta_V \Omega_{\Gamma}}\bigg)^{\frac{1}{2}} \big( \beta_\Gamma e^{-i\Omega_\Gamma t} +\beta^\dagger_{\Gamma} e^{i\Omega_\Gamma t} \big).
\end{equation}

Using the symmetry point group analysis, we identify the Raman-active phonon modes by considering only the vibrations of the Iridium atoms:

\begin{equation}
\Gamma^{\text R}_\text{Ir}=A_g\oplus B_{1g}\oplus2B_{2g}\oplus2B_{3g}.
\end{equation}

Details  of the analysis are provided in Appendix <ref>.

As mentioned in the introduction, the Fano lineshape has been observed experimentally in $\beta$-Li$_2$IrO$_3$ for the low-energy $A_g$  phonon mode centered around 24 meV <cit.>.
This $A_g$ phonon involves out-of-phase vibrations of Iridium atoms in the $c$-direction and stretching vibrations of the Oxygen octahedral cage <cit.>. We visualize this mode in Fig. <ref>, where we focus on the vibrations of the Iridium atoms and predominantly ignore the vibrations of the oxygen atoms, as these mainly contribute to non-Kitaev interactions. The decomposition of the phonon displacement field $u_{A_g}$ is detailed in Appendix <ref>.

The bare phonon propagator in imaginary time $\tau=it$ for a given irrep $\Gamma$ is defined as $\mathcal{D}_\Gamma^{(0)}(\tau) = - \braket{T_\tau u_\Gamma(\tau)u_{\Gamma}(0)}$, where $T_τ$ is the imaginary time ordering operator.  Performing the Fourier transform to Matsubara frequencies,
$𝒟_Γ^(0)(iΩ_m)=∫_0^βdτ e^iΩ_m τ𝒟_Γ^(0)(τ)$, we obtain
\begin{equation}
\mathcal{D}^{(0)}_\Gamma(i\Omega_m) = -\bigg( \frac{\hbar}{\rho\delta_V} \bigg) \frac{1}{(i\Omega_m)^2-\Omega_\Gamma^2}.
\end{equation}

\begin{figure}[H]
\centering
\includegraphics[width=0.75\linewidth]{fig2.pdf}
\caption{Visualization of the $24$ meV  $A_{g}$ optical phonon mode. We only consider the vibration of 4 Iridium atoms (shown in red) to preserve the exact solvability of the model and capture the underlying dynamics with the pure Kitaev model. The group theoretical analysis used to obtain this mode is presented in Appendix \ref{subsec:eigenmode}.}
\label{fig:AgPhonon}
\end{figure}

\subsection{Magnetoelastic coupling}\label{sec:sp-phCoupling}

Lattice vibrations can lead to changes in bond lengths and bond angles, which in turn affect the Kitaev interaction strength $J_K$.
Although these vibrations can potentially introduce new couplings within the system, such as Heisenberg, symmetric $Γ$, and antisymmetric Dzyaloshinskii-Moriya interactions, we will neglect these effects here, assuming they are small.
Our focus will be on capturing the phonon dynamics while preserving the solvability of the Kitaev spin model.

From a microscopic perspective, the coupling between the Kitaev QSL and phonon modes can be written as~\cite{Ye2020}: \begin{equation}
\mathcal{H}_{\text{s-ph}}=-\lambda\sum_{\mathbf{r},t}\mathbf{M}_t \cdot [\mathbf{u(r)-u(r+M}_t)]\sigma_\mathbf{r}^{\alpha_t}\sigma_{\mathbf{r+M}_t}^{\alpha_t}, %\lambda\sum_{\mathbf{r},\alpha}\mathbf{M}_{\alpha}\cdot [\mathbf{u(r)-u(r+M_\alpha)}]\sigma_\mathbf{r}^{(\alpha)}\sigma_\mathbf{r+M_\alpha}^{(\alpha)}
\end{equation}
where index $t$  runs over the nearest neighbor bonds  shown in Fig.~\ref{fig:3dnn} and $u(r)$ is the lattice displacement vector at the lattice point $r$.
%$\mathbf{M_{x}}$,  $\mathbf{M_{x'}}$, $\mathbf{M_{y}}$, $\mathbf{M_{y'}}$, and  $\mathbf{M_{z}}$ denote the relevant lattice vectors shown in Fig.~\ref{fig:3dnn}.
By rewriting the spin bilinears in terms of Majorana fermions, we obtain a coupling between phonons and the fractionalized excitations of the Kitaev QSL.
In the Raman response, the fractionalized excitations of the QSL form a continuum, which, through the discussed magneto-elastic coupling, can result in a characteristic Fano lineshape of the optical phonons.
As previously discussed, the Fano lineshape for the 24 meV $A_g$ phonon in $β$-Li$_2$IrO$_3$ has been observed, and understanding its spectral features is the focus of this study.

We can also express
the spin-phonon coupled Hamiltonian using the symmetry point group $D_2h$ in the $A_g$ irrep as
\begin{equation}
\mathcal{H}_{\text{s-ph}}^{A_g}= u_{A_g}\big[2\lambda_{A_g}^{zz}\Sigma_{A_g}^{zz}-\lambda_{A_g}^{xx}\Sigma_{A_g}^{xx}-\lambda_{A_g}^{yy}\Sigma_{A_g}^{yy}\big],
\end{equation}
\begin{equation}
\begin{aligned}
\Sigma_{A_g}^{xx}&=\sum_{\mathbf{r}}\bigg(\sigma_\mathbf{r}^x\sigma_\mathbf{r+M_x}^x+\sigma_\mathbf{r}^x\sigma_\mathbf{r+M_{x\prime}}^x\bigg)\\
\Sigma_{A_g}^{yy}&=\sum_{\mathbf{r}}\bigg(\sigma_\mathbf{r}^y\sigma_\mathbf{r+M_y}^y+\sigma_\mathbf{r}^y\sigma_\mathbf{r+M_{y\prime}}^y\bigg)\\
\Sigma_{A_g}^{zz}&=\sum_{\mathbf{r}}2\,\sigma_\mathbf{r}^z\sigma_\mathbf{r+M_z}^z.
\end{aligned}
\end{equation}
If we assume isotropic spin-phonon coupling strengths, $λ_A_g^xx∼λ_A_g^yy∼λ_A_g^zz=λ_A_g$,  the spin-phonon coupling Hamiltonian can be   simplified as
\begin{equation}\label{HsphA1g}
\mathcal{H}_{\text{s-ph}}^{A_g}= \lambda_{A_g}u_{A_g}\big[2\Sigma_{A_g}^{zz}-\Sigma_{A_g}^{xx}-\Sigma_{A_g}^{yy}\big].
\end{equation}
To rewrite  the spin-phonon Hamiltonian in the Bogoliubov basis,
we first obtain the spin-phonon coupling vertex $Λ_𝐤$ in the Majorana fermion representation and then perform the transformation to the Bogoliubov basis.
In the Majorana basis, the  spin-phonon coupling matrix in
the momentum space reads:
\begin{equation}
\label{eqn:lambdaK}
\Lambda_\mathbf{k}= \begin{pmatrix}
0 & 0 & \Tilde{B}_{\mathbf{k}} & \Tilde{A}_{\mathbf{k}} \\
0 & 0 & \Tilde{A}_{\mathbf{k}} & \Tilde{C}_{\mathbf{k}} \\
\Tilde{B}^*_{\mathbf{k}} & \Tilde{A}^*_{\mathbf{k}} & 0 & 0 \\
\Tilde{A}^*_{\mathbf{k}} & \Tilde{C}^*_{\mathbf{k}} & 0 & 0
\end{pmatrix},
\end{equation}
where %$N$ is the number of unit cells and
$A_𝐤=-2i J_K$, $B_𝐤=i J_K( e^i𝐤·𝐚_2+ e^i𝐤·𝐚_1)$,  and $C_𝐤=iJ_K(1+ e^-i𝐤·𝐚_3)$.
Then we
use the unitary  transformation matrix  $U_k$ from Eq.~\eqref{eqn:c2a}, and get
\begin{equation}\Tilde{\Lambda}_{\bf k}=U^\dagger_{\bf k}\Lambda_{\bf k} U_{\bf k}.
\end{equation}
Now, the spin-phonon Hamiltonian can be written as
\begin{equation}\label{eqn:sp-phHamAg}
\mathcal{H}_{\text{s-ph}}^{A_g}=\lambda_{A_g}\sum_\mathbf{k}A_\mathbf{k}^\dagger\Tilde{\Lambda}_\mathbf{k}A_\mathbf{k}.
\end{equation}

%{\cbl commutation of spin-phonon coupling matrix with hamiltonian matrix }
The spin-phonon Hamiltonian (\ref{eqn:sp-phHamAg}) does not commute with the spin Hamiltonian (\ref{eqn:spinhambogoliu}), resulting in a non-zero renormalization of the phonon propagator, as discussed in the subsequent sections.
This is a consequence of the independent basis functions  in the $A_g$ irrep of the point group $D_2h$ .

\subsection{Phonon Polarization Bubble}\label{sec:polbub}

To study the effects of spin-lattice coupling, we evaluate one-loop corrections to the phonon self-energy, represented by the phonon polarization bubble.
In the finite temperature formalism, this is determined in imaginary time as:
\begin{equation}\label{eqn:polbubexpr}
\Pi(\tau)=\frac{\lambda_{A_g}^2}{N}\sum_\mathbf{k}\braket{(T_\tau A_\mathbf{k}^\dagger\Tilde{\Lambda}_\mathbf{k} A_\mathbf{k})(\tau)(T_\tau A_\mathbf{k}^\dagger\Tilde{\Lambda}_\mathbf{k} A_\mathbf{k})(0)},
\end{equation}
where $N$ denotes the number of units  cells.
Details of the derivation of phonon polarization bubble and relevant notation can be found in Appendix~\ref{Appendix:polbub}.
The final expression obtained in Matsubara space is presented below:
\begin{subequations}\label{eqn:polbubres}
\begin{align}
\Pi^{\text{pp}}(i\Omega_m)=\frac{\lambda_{A_g}^2}{N}\sum_\mathbf{k}&\bigg(P_{aa}^{gg}\Tilde{\Lambda}_\mathbf{k}^{23}\Tilde{\Lambda}_\mathbf{k}^{32} + P_{bb}^{gg}\Tilde{\Lambda}_\mathbf{k}^{14}\Tilde{\Lambda}_\mathbf{k}^{41}
\\&+P_{ba}^{gg}\Tilde{\Lambda}_\mathbf{k}^{13}\Tilde{\Lambda}_\mathbf{k}^{31} + P_{ab}^{gg}\Tilde{\Lambda}_\mathbf{k}^{24}\Tilde{\Lambda}_\mathbf{k}^{42}
\bigg)\nonumber,\\
\begin{split}\label{eqn:polbubresph}
\Pi^{\text{ph}}(i\Omega_m)=\frac{\lambda_{A_g}^2}{N}\sum_\mathbf{k}  &P_{ba}^{g\Bar{g}}\bigg(\Tilde{\Lambda}_\mathbf{k}^{34}\Tilde{\Lambda}_\mathbf{k}^{43} +\Tilde{\Lambda}_\mathbf{k}^{12}\Tilde{\Lambda}_\mathbf{k}^{21}
\bigg),
\end{split}
\end{align}
\end{subequations}
where pp- and ph- refers, correspondingly,  to
the particle-particle and the particle-hole  channels
of phonon scattering from fractionalized fermionic excitations of the Kitaev spin liquid. We introduce $P_γγ'^gg(iΩ_m)$ and $P_γγ'^gg̅(iΩ_m)$ as a shorthand for the convolution of Matsubara Green's functions explicitly shown in Eq.~\eqref{eqn:matsum}. For brevity, the Matsubara frequency dependence has been suppressed in the above equations. It is also worth noting that intraband processes
in the ph-channel only contribute to Rayleigh scattering as can be seen from the lack of $P_aa^gg̅$ and $P_bb^gg̅$ type terms in Eq.~\eqref{eqn:polbubresph}.
The sum over momentum space is performed numerically, and the real and imaginary parts of the polarization bubble are plotted in

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{fig3.pdf}
\caption{Plot of the real and imaginary parts of the polarization bubble $\Pi(\Omega)$, measured in units of  $\lambda_{A_g}^2$, as a function of frequency, $\Omega$ measured in units of $J_K$, computed  at  $\text{T}/J_K=0.1$. The dashed line is marking  the frequency of the optical phonon.}
\label{fig:piplt}
\end{figure}

The polarization bubble renormalizes the phonon propagator via the Dyson equation as $𝒟=[(𝒟^(0))^-1-Π]^-1$, where $Π=Π'+i  Π”$ is expressed in terms of its real and imaginary parts. In particular,  the phonons acquire a finite life-time given by
the imaginary part of the renormalized phonon propagator:
\begin{equation}\label{eq:broadened_phonon}
\mathrm{Im\, \mathcal{D}}(\Omega)\!=\!\dfrac{\lambda_{A_g}^2 \,\Pi''(\Omega_{A_g})}{\left[%\frac{\delta_V\rho}{\hbar}
\Omega_{A_g}^2\!-\!\Omega^2\!-\!\lambda_{A_g}^2 \,\Pi'(\Omega_{A_g})\right]^2\!+\!\lambda_{A_g}^4 \,\Pi''(\Omega_{A_g})^2}.
\end{equation}
Furthermore, as will be seen later, the phonon Raman peak parameters, such as the width, center position, and asymmetry factor, are directly related to the temperature dependence of $Π'$ and $Π”$ (Fig. \ref{fig:pitmpplt}).

\section{Methods}\label{sec:method}

In this section, we present the diagrammatic perturbation theory framework used to evaluate Raman intensity, as developed in Ref.~\cite{feng2022}.
First, we construct the vertices relevant to the scattering processes that obey the kinematic constraints. Using these vertices, we then construct the intensity diagrams and evaluate the Raman response.

\subsection{ Magnetic  Raman  operator}\label{sec:spin-raman}

The magnetic Raman  operator describes the interaction between photons and spin degrees of freedom. According to the Loudon-Fleury theory \cite{fleury1968scattering}, it can be  expressed as
\begin{equation}\label{RamanLF}
\mathcal{R}_{\text{em-s}}=\sum_{\mu,\mu'}\varepsilon^\mu_{\text{in}}\mathcal{R}^{\mu \mu'}_{\text{em-s}}\varepsilon^{\mu'}_{\text{out}},
\end{equation}
\begin{equation}
\mathcal{R}^{\mu \mu'}_{\text{em-s}}=\nu\sum_{\mathbf{r},t}\mathbf{M}_t^\mu \mathbf{M}_t^{\mu'} \sigma_\mathbf{r}^{\alpha_t}\sigma_{\mathbf{r+M}_t}^{\alpha_t}.
\end{equation}
Here $ε_in$ and $ε_out$ are the  polarization vectors of the incoming and outgoing light and $ν$ is the photon-spin coupling constant.
can also be decomposed into the irreps of the point group as follows:
\begin{equation}\mathcal{R}^{\mu \mu'}_{\text{em-s}}=\sum_{\Gamma} \alpha_\Gamma R_\Gamma^{\mu\mu'}\Sigma_{\Gamma},
\end{equation}
where the coefficient $α_Γ=1/2 Tr[R_Γ^T ·ℛ_em-s]$ is obtained using orthogonality of irreps, and $R_Γ$ is Raman tensor of irrep $Γ$. Here, we are interested in the magnetic Raman  operator in the  $A_g$ irrep of $D_2h$. Due to the independent nature of the basis functions in this irrep, the Raman operator can be written as  $ℛ^A_g_em-s=ℛ^aa_em-s+ℛ^bb_em-s+ℛ^cc_em-s$.

\subsection{ Phonon Raman operator}\label{sec:phonon-raman}
The phonon Raman  operator is determined by the interaction between the phonons and the polarization vectors of the incoming and outgoing light.
The phonon Raman operator in the $A_g$ irrep is   given by
\begin{equation}\mathcal{R}^{A_g}_{\text{em-ph}}=\mu_{A_g} \mathcal{R}_{A_g}u_{A_g},
\end{equation}
where the Raman tensor is given by \cite{bilbaoSAM}:
\[
a & 0 & 0\\
0 & b & 0\\
0 & 0 & c
\end{pmatrix}
\]
For  our four-sublattice structure, the Raman tensor is simply  given by $ℛ_A_g=ℐ_4⊗R_A_g$, where $ℐ_4$ is the $4×4$ identity matrix.  In the following, we set  $a=b=0$, $c=1$ to consider $cc$ polarization as  experimentally reported in \cite{glamazda2016expt}.

\subsection{ Raman intensity: Fano shape of the  24 meV  $A_g$ phonon }\label{sec:fano}

Knowing the Raman operator, we can compute the Raman intensity as $I(τ)=-
⟨T_τ ℛ(τ)ℛ(0)⟩,$
where  $⟨⋯⟩$ is the thermal  average and $T_τ$ is the imaginary time ordering operator.
The effect  of  the spin-phonon coupling  modifies the Raman response, and it  can be computed systematically
by performing the $S-$matrix expansion and expressing the coupling in the interaction picture. For finite temperatures, the intensity is given by:
\begin{equation}
I(\tau)=-\braket{T_\tau\,\mathcal{R}(\tau)\mathcal{R}(0)e^{-\int_0^\beta d\tau'\mathcal{H}_{\text{s-ph}}(\tau')}},
\end{equation}
where $ℛ=ℛ_em-s+ℛ_em-ph$ and only distinct connected diagrams are summed over.  The perturbative expansion yields the following expression:
\begin{equation}
I(\tau)=\sum_{k=0}^\infty \, (-1)^k\,\prod_i^k\int_0^\beta d\tau_i \, \braket{T_\tau \mathcal{R}(\tau)\mathcal{R}(0)\prod_i^k\mathcal{H}^{A_g}_{\text{s-ph}}(\tau_i)}.
\end{equation}
Performing the Fourier transform to the frequency domain,
we get:
\begin{equation}
I(i\Omega_m) =  -\int_0^\beta  d\tau e^{i\Omega_m \tau} \,I(\tau),%\langle T_\tau  \mathcal{R}({\tau})\ \mathcal{R}(0)\rangle,
\end{equation}
where $β=1/T$
(with $k_B=1$ being the Boltzmann constant).

At zeroth order, $k=0$, we get the pure contribution from the magnetic and phonon Raman intensities, $I^(0)_A_g=I^A_g_em-ph+I^A_g_em-s$.
%Part of this intensity arises from purely magnetic and phonon responses.
The component of the magnetic part of the intensity is given by:
\begin{equation}
I_{\text{em-s}}^{A_g}(i\Omega_m)=  -\int_0^\beta  d\tau e^{i\Omega_m \tau}\braket{T_{\tau}\mathcal{R}^{A_g}_{\text{em-s}}(\tau)\mathcal{R}^{A_g}_{\text{em-s}}(0)}.
\end{equation}
The explicit form of the Raman matrix and the intensity can be found in Appendix \ref{app:mag-raman}.
%{\cred NP: next sentence is not clear here, cc should be explained or removed}\\
The symmetry arguments for the  Loudon-Fleury form of the Raman operator give $I^cc=9I^aa$ and $I^bb=4I^aa$~\cite{perreault2015theory}. Thus,   we have $I_em-s^A_g ∝I^cc$ corresponding to the polarization geometry in the experimental setup for the Raman scattering in  $β$-Li$_2$IrO$_3$ \cite{glamazda2016expt}.
The phonon contribution to the Raman intensity is calculated as:
\begin{equation}
I_{\text{em-ph}}^{A_g}(i\Omega_m)=  -\int_0^\beta  d\tau e^{i\Omega_m \tau}\braket{\mathcal{R}^{A_g}_{\text{em-ph}}(\tau)\mathcal{R}^{A_g}_{\text{em-ph}}(0)},
\end{equation}
which yields a sharply peaked delta function at the phonon frequency $Ω_A_g$.

At first order, $k=1$, the  contributions to the Fano asymmetry arise. The intensity can be written as:
\begin{equation}
I_{A_g}^{(1)}(\tau) = \int_0^\beta d\tau_1 \,\braket{T_\tau\mathcal{R}(\tau)\mathcal{R}(0)\mathcal{H}^{A_g}_{\text{s-ph}}(\tau_1)}.
\end{equation}
It can be further decomposed as
\begin{equation}\label{eqn:fanoInt}
I_{A_g}^{(1)}(\tau) =I_L^{(1)}(\tau)+I_R^{(1)}(\tau),
\end{equation}
where $I_L^(1)(τ)=∫_0^β  dτ_1  ⟨T_τℛ^A_g_em-ph(τ)ℛ^A_g_em-s(0)ℋ^A_g_s-ph(τ_1)|$⟩, and $I_R^{(1)}(\tau)=\int_0^\beta\, d\tau_1 \,\braket{T_\tau\mathcal{R}^{A_g}_{\text{em-s}}(\tau)\mathcal{R}^{A_g}_{\text{em-ph}}(0)\mathcal{H}^{A_g}_{\text{s-ph}}(\tau_1)}$. These contributions are diagrammatically represented in  Eq.~\eqref{eqn:FanoDiags} where panel (a) is $I_L^(1)(τ)$ and panel (b) is $I_R^(1)(τ)$:

\iffalse
\begin{figure}[H]
\includegraphics[width=\linewidth]{images/FanoDiags.pdf}
\caption{Diagrams contributing to the Fano asymmetry in the Raman intensity corresponding to the terms in Eq.\eqref{eqn:fanoInt}. Panel (a) is $I_R^{(1)}(\tau)$ and panel (b) is $I_L^{(1)}(\tau)$}
\label{fig:FanoDiags}
\end{figure}
\fi
\begin{equation}\label{eqn:FanoDiags}
\begin{gathered}
\includegraphics[width=0.65\linewidth]{eqn_fig.pdf}
\end{gathered}
\end{equation}
Here, the wavy lines represent photon propagators, the dashed lines denote phonon propagators renormalized by the spin-phonon interaction, and the solid lines indicate Majorana fermion propagators. The details for the evaluation of both the left and the right diagrams are outlined in Appendix \ref{App:Fano}.
It is also important to note that all fermionic lines in the intensity diagrams originate from Majorana fermions, as we do not account for the $ℤ_2$ fluxes in our consideration of various scattering processes.  This is because $ℤ_2$ fluxes are gapped, and there is a finite-temperature phase transition separating the zero-flux phase from the finite-temperature phase with flux excitations \cite{eschmann2020thermodynamic}.

Finally, we perform a Fourier transform to the Matsubara frequency, followed by an analytical continuation to real frequency, where $i Ω_m →Ω+ iδ$. Here, $Ω= Ω_in - Ω_out$ represents the frequency shift, defined as the difference between the incoming and outgoing photon frequencies. The total intensity is then given by:
\begin{equation}
I_{A_g} (\Omega)=I_{A_g}^{(0)} (\Omega)+I_{A_g}^{(1)}(\Omega).
\end{equation}

\section{Results}\label{sec:results}
To compare our theoretical results with experimental data, we need to fix the model parameters that appear in the total intensity. We detail this procedure in the next subsection, followed by a presentation of our findings.

\subsection{Fitting the model parameters}\label{sec:fitparams}

There are  free  coupling parameters $(μ_A_g,ν,λ_A_g)$ that need to be fixed before performing further quantitative computations.
First, we note that $μ_A_g$ and $ν$ control the overall intensity of the phonon Raman peak and the magnetic continuum, respectively. Since the overall magnitude of the Raman spectrum is free to rescale, increasing or reducing $μ_A_g$ and $ν$ uniformly by the same factor will retain the shape of the calculated spectrum, differing only in overall scale.
In contrast, $λ_A_g$, the spin-phonon coupling strength,   cannot be freely rescaled as it controls the phonon peak widths as well as its Fano asymmetry. Thus, fixing $λ_A_g$ with respect to $μ_A_g$ and $ν$   is crucial for accurately determining the detailed features of the overall Raman spectrum.

We fit these model parameters at a low temperature using the experimental data in \cite{glamazda2016expt}. We estimate $ν$ from the ratio of the peak to continuum reported in this study, giving
$ν^2/μ^2∼1/6$, which suggests $ν∼0.41 J_K$. To set $λ_A_g$, we consider the imaginary part of the renormalized phonon propagator, expected to be a symmetric Lorentzian depending only on $λ_A_g$. By fitting it to the experimental peak with the asymmetry manually removed, we estimate $λ_A_g=0.49 J_K$.
Then, we fix the overall scale by setting
$μ_A_g=1,$ and further fine-tune these model parameters by comparing the sum of all relevant diagrams to the experimental phonon Raman peak. A comparison between the experimental data and our calculation is shown in Fig.~\ref{fig:modelparam} with the fitted parameter values $(μ_A_g,ν,λ_A_g)=(1.0,-0.37,0.29)$ in units of the Kitaev strength $J_K$.

\begin{figure}[t]
\centering
\includegraphics[width=0.95\linewidth]{fig4.pdf}
\caption{Plot of the experimental phonon Raman peak (dashed blue) compared with the theoretically evaluated Raman intensity (solid red) for the model parameters $\mu_{A_g}=1.0, \nu=-0.37$ and $\lambda_{A_g}=0.29$. All coupling  are in the units of $J_K$.}
\label{fig:modelparam}
\end{figure}

The Fano lineshape of the Raman intensity at any temperature can be phenomenologically captured by the asymmetric Lorentzian function:
\begin{equation}\label{eqn:asymLor}
\end{equation}
Here, the parameters $I_0,q,ω_0$, and $γ$ control the overall intensity, asymmetry, peak position, and peak width respectively. We present the temperature evolution of the total Raman intensity in terms of the asymmetry parameter, $1/|q|$, the peak position, $ω_0$, and the peak width, $γ$, in the following section.

\subsection{Temperature Evolution of the Fano Peak Parameters}\label{sec:tempevol}
Once the model parameters $μ_A_g$, $ν$, $λ_A_g$ are fixed we calculate the Raman intensity,  $I(Ω)$, as a function of temperature. We then fit our calculated intensity to the asymmetric Lorentzian Eq.~\eqref{eqn:asymLor} to extract the peak parameters at each temperature. The temperature evolution of the calculated Fano peak parameters along with the experimental data \footnote{We divide by 2 the linewidth shown in temperature evolution data presented in Ref. \cite{glamazda2016expt} to obtain the correct FWHM for comparison. $\Gamma$ represents the FWHM where $\gamma$ is the HWHM such that $\Gamma=2\gamma$.} is shown in Fig.~\ref{fig:peakparamtemp}. All momentum-space integrals are performed numerically. The momentum-space is discretized into $100×100×100$ points. The frequency-space is discretized into $201$ points.

\begin{figure}[t]
\centering
\includegraphics[width=0.95\linewidth]{fig5.pdf}
\caption{Theoretically computed temperature evolution of the Fano peak parameters plotted with the experimental data (shown by black circles) \cite{glamazda2016expt}. Blue ($T_L=27\,K$) and red ($T_H=209\,K$)  dashed lines mark temperature scales characteristic of the Kitaev model \cite{do2017majorana}. The spin-phonon contribution is shown in solid colors, with the grayed-out part past $T_H=209\,K$ indicating the loss of validity of the spin-phonon model as spins are no longer fractionalized.  The phonon-phonon contribution, shown in open circles, accounts for the anharmonic decay of optical phonons as described by
Eq.~\eqref{eqn:phononAnhar}. Phenomenological anharmonic model constants are $(A,B,C,D) = (-7.99,\,-0.0198,\,19.7,\,0.0218) \mu \mathrm{eV}$. We use Kitaev strength of $J_K=18$ meV for unit conversions.}
\label{fig:peakparamtemp}
\end{figure}

A note is in order here. Phonon scattering from fractionalized excitations occurs at relatively low temperatures, with the range determined by the magnitude of the Kitaev interaction. At higher temperatures, phonon-phonon scattering processes become dominant. Although we do not consider these explicitly, we  account for phonon-phonon scattering processes as described by the phenomenological anhamonicity model \cite{balkanski1983anharmonic}:
\begin{equation}\label{eqn:phononAnhar}
\begin{aligned}
\omega(T) = \omega_0 + C\left[1 + \frac{2}{e^{x} - 1}\right] + D\left[1 + \frac{3}{e^{y} - 1} + \frac{3}{(e^{y} - 1)^2}\right]\\
\Gamma(T) = \Gamma_0 + A\left[1 + \frac{2}{e^{x} - 1}\right] + B\left[1 + \frac{3}{e^{y} - 1} + \frac{3}{(e^{y} - 1)^2}\right]
\end{aligned}
\end{equation}
where $x=ħω_0/2k_B T$,  $y=ħω_0/3k_B T$ and $ω_0$ and $Γ_0$ are the frequency and the linewidth of an optical mode at very low temperature, respectively. Supplementing our spin-phonon coupling model with the phonon-phonon decay at higher temperatures provides a good explanation for the experimental data.

Fig.~\ref{fig:peakparamtemp} (a) shows the temperature evolution of the Fano asymmetry, with purple circles for calculated values and black circles for experimental values. The Fano asymmetry decreases as temperature rises due to reduced pp-scattering from occupied low-energy fermionic states.

Fig.~\ref{fig:peakparamtemp} (b) shows the temperature evolution of $Γ$ (FWHM) of the total Raman intensity. Experimental values are shown as filled black circles, spin-phonon model calculations as filled blue circles, and anharmonic phonon-phonon model calculations as open blue circles. The decrease in $Γ$ with increasing temperature is due to reduced pp-scattering, leading to less broadening.

Fig.~\ref{fig:peakparamtemp} (c) shows the temperature evolution of the peak position $ω_0$ of the total Raman intensity. Experimental values are shown as filled black circles, spin-phonon model calculations as filled red circles and anharmonic phonon-phonon model calculations as open red circles. The increase in $ω_0$ with increasing temperature can be understood from the expression for broadened phonon peak given in Eq.~\eqref{eq:broadened_phonon}. Since the real part of the phonon polarization bubble $Π'$ decreases with temperature the peak of the phonon Raman intensity shifts to higher values as can be seen from the denominator of Eq.~\eqref{eq:broadened_phonon} and Fig. \ref{fig:pitmpplt}. %\textcolor{red}{[Maybe we need the temperature evolution of Im, Re parts of $\Pi$ in Appendix.]}

As expected, the agreement between experimental and calculated values from the spin-phonon coupled model in  Fig.~\ref{fig:peakparamtemp} (a)-(c)  is better for temperatures below $T ∼150$ K, marked by the dashed black grid-line.  %{\cred   NP: I do not think that we can be that specific about 156 K. Shall we just say above 150K?}
Below this temperature, the Majorana excitations of the Kitaev QSL  almost solely determine the shape of the 24 meV phonon mode through magnetoelastic coupling, resulting in a better fit with the experimental data.  However, above $T ∼150$ K, anharmonic effects become dominant and show better agreement with the experimental data, as seen in Fig.~\ref{fig:peakparamtemp} (b), (c). Additionally, flux proliferation with increasing temperature would contribute to a further decrease in the asymmetry characteristic of the Fano lineshape.

\subsection{$J_K$ dependence of Fano asymmetry}\label{sec:jkdep}

We now discuss the dependence of the Fano asymmetry, characterized by $1/|q|$, on the Kitaev coupling strength $J_K$, as shown in Fig.~\ref{fig:jkvar}.
The asymmetry of the Fano lineshape is determined  by the relative position of the phonon peak and the Majorana fermion continuum.
While the spread of the continuum scales with $J_K$, the position of phonon Raman peak depends only on the
elastic constants. Therefore, a change in $J_K$ modifies the relative position of the phonon peak to the Majorana continuum, leading to a change in $1/|q|$.

In Fig.~\ref{fig:jkvar}, we observe this  $J_K$ dependence of the Fano asymmetry through the temperature evolution of $1/|q|$ for three values of
$J_K=12, 18,$ and $24$ meV, shown in green, purple, and blue, respectively,with experimental data represented by filled black circles. Note that while all three curves follow the same trend of monotonically decreasing with temperature, the Fano asymmetry is suppressed more for larger values of
$J_K$. This trend can be understood from the aforementioned points. As $J_K$ increases the weight of the Majorana continuum shifts to the right of the phonon peak leading to smaller asymmetry.

This dependence of the Fano asymmetry on $J_K$  can help in estimating the strength of the Kitaev coupling in
$β$-$Li_2IrO_3$. A simple comparison of the obtained curves with the experimental findings in Fig.~\ref{fig:jkvar}  suggests that 18 meV is a very plausible Kitaev coupling strength, further confirming our previous estimates from other dynamical responses \cite{ducatman2018magnetic,Rousochatzakis2018,Li2020,Li2020b,Yang2022,Halloran2022,Ruiz2021}. %\textcolor{teal}{It also shows that $J_K$ is bounded below 24 meV.}

\begin{figure}[t]
\centering
\includegraphics[width=0.95\linewidth]{fig6.pdf}
\caption{ The Fano parameter $1/|q|$ as a function of temperature, shown for three Kitaev couplings $J_K$ within the range of proposed Kitaev strengths for $\beta$-Li$_2$IrO$_3$. The Fano asymmetry is significantly suppressed at higher Kitaev strengths, setting a limit for the possible values of $J_K$.
\label{fig:jkvar}
\end{figure}

\section{Conclusion}\label{sec:conclucsions}

In this work, we presented a theoretical investigation aimed at explaining the temperature dependence of the Fano lineshape of the 24 meV phonon peak observed by Raman spectroscopy in the three-dimensional Kitaev candidate material  $β$-$Li_2IrO_3$ \cite{glamazda2016expt}. We achieved this by considering a pure Kitaev model on a 3D hyperhoneycomb lattice, coupled with low-energy optical phonons.  In this model, spins fractionalize into Majorana fermions and $Z_2$ fluxes,
However, at sufficiently low temperatures, below the flux phase transition \cite{eschmann2020thermodynamic}, the spin-phonon  interaction is between the phonons and the Majorana fermions. Therefore, we assumed that the 24 meV phonon is predominantly coupled to Majorana fermions at low temperatures. As the temperature increases and fluxes become thermally excited, the shape of the phonon peak is already mainly determined by anharmonicity effects. We expect fluxes to only add small corrections to this behavior, which we choose to disregard in the current work due to the complexity they introduce in the numerics.

We find that our spin-phonon coupling model closely captures the temperature evolution of the Fano asymmetry and other peak parameters up to the major fermionic excitation temperature of approximately 150 K. Beyond this temperature, Majorana fermions are no longer good degrees of freedom, and the scattering processes become dominated by phonon-phonon anharmonicities.  By combining these two mechanisms, we obtained  very good agreement with the experimental results \cite{glamazda2016expt}.   We also used our formalism to refine the range of the Kitaev coupling strength in $β$-Li$_2$IrO$_3$. By analyzing the Fano asymmetry parameter $1/|q|$, we estimated that the Kitaev coupling strength is bounded between 18 meV and 24 meV.

To conclude, our study of the Fano lineshape in Raman spectroscopy provides valuable insights into the temperature dependence of the 24 meV phonon peak in $β$-Li$_2$IrO$_3$ and demonstrates the potential of phonon Raman spectroscopy as a probe for other quantum spin liquids. Our theoretical approach to analyzing the Fano asymmetry and its dependence on coupling strengths, temperature, and other external parameters can be easily  extended to various spin liquid candidates. This makes it a powerful tool for extracting information about quantum spin liquid phases from Raman studies of the phonon Fano lineshape.

\section*{Acknowledgements}
The authors thank  P. Peter Stavropoulos, Vitor Dantas, Kexin Feng, Yang Yang, David Mayrhofer and Andrey Chubukov for valuable discussions. The work is supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0018056.   We acknowledge the support from NSF DMR-2310318 and  the support of the Minnesota Supercomputing Institute (MSI) at the University of Minnesota. N.B.P. also acknowledges the hospitality and partial support  of the Technical University of Munich – Institute for Advanced Study and  the Alexander von Humboldt Foundation.\\

\appendix

\section{Fractionalization  in the Kitaev model on the hyperhoneycomb lattice }\label{appendix:geometry}

In this Appendix, we review the solution of the Kitaev model on the hyperhoneycomb lattice, earlier presented in Refs.~\cite{OBrien2016,perreault2015theory,halasz2017RIXS,feng2022sound}.

The primitive unit cell of the hyperhoneycomb lattice (see Fig.~\ref{fig:3dnn}) comprises of four sublattices indexed by $δ=1,2,3,4$. Their coordinates in the orthorhombic basis are:
\begin{eqnarray}
\delta_2=(0,0,\frac{c}{6}),\\
\delta_4=(\frac{a}{4},\frac{b}{4},\frac{5c}{12}).\nonumber
\end{eqnarray}
The orthorhombic unit cell contains
four  primitive unit cells. The lattice vectors of the orthorhombic unit cell are
\begin{eqnarray}
\mathbf{a}=(a,0,0),\,
\mathbf{b}=(0,b,0),\,\mathbf{c}=(0,0,c),
\end{eqnarray}
such that $a:b:c=1:√(2):3$. The Cartesian axes are related to the orthorhombic axes as
\begin{eqnarray}
\mathbf{x}=(\mathbf{a}+\mathbf{c})/\sqrt{2},\,\mathbf{y}=(\mathbf{c}-\mathbf{a})/\sqrt{2},\, \mathbf{z}=-\mathbf{b}.
\end{eqnarray}
The primitive unit cell lattice vectors are then given  by
\begin{eqnarray}
\mathbf{a_1}=(\frac{a}{2},0,\frac{c}{2}),\mathbf{a_2}=(0,\frac{b}{2},\frac{c}{2}),\mathbf{a_3}=(\frac{a}{2},\frac{b}{2},0),
\end{eqnarray}
so that
the reciprocal lattice vectors are
\begin{eqnarray}
\end{eqnarray}
In this convention, the spin-Hamiltonian can be represented as a matrix:
\begin{widetext}
\begin{eqnarray}
\mathcal{H}_{\text{s}}=\sum_\mathbf{k} C^\dagger_\mathbf{k}H_\mathbf{k}^{s}C_\mathbf{k}=\sum_{\mathbf{k}}
\begin{pmatrix}
c_{\mathbf{-k},4} & c_{\mathbf{-k},2} & c_{\mathbf{-k},1} & c_{\mathbf{-k},3}
\end{pmatrix} \begin{pmatrix}
0 & 0 & B_{\mathbf{k}} & A_{\mathbf{k}} \\
0 & 0 & A_{\mathbf{k}} & C_{\mathbf{k}} \\
B^*_{\mathbf{k}} & A^*_{\mathbf{k}} & 0 & 0 \\
A^*_{\mathbf{k}} & C^*_{\mathbf{k}} & 0 & 0
\end{pmatrix} \begin{pmatrix}
\end{pmatrix},\label{H_matrix}
\end{eqnarray}
where $A_{\mathbf{k}}=-iJ_K$, $B_{\mathbf{k}}=-i J_K(\,e^{i\mathbf{k}\cdot\mathbf{a_2}}+\,e^{i\mathbf{k}\cdot\mathbf{a_1}})$, and $C_{\mathbf{k}}=-i J_K(1+\,e^{-i\mathbf{k}\cdot\mathbf{a_3}})$.
\end{widetext}
The quadratic fermionic Hamiltonian (\ref{H_matrix})
can be diagonalized via a standard procedure. Since the hyperhoneycomb lattice has four sites per
unit cell, the resulting band structure has four fermion
bands shown in Fig.~\ref{fig:ebandstruct}.  Note that  due to the particle-hole redundancy of the   Kitaev model in terms of Majoranas,
the physical  excited states are described  by
the two positive bands only. Furthermore, the lowest positive band displays a nodal line structure \cite{OBrien2016,ybk2015prl}.

\begin{figure}[h]
\centering
\includegraphics[width=0.95\linewidth]{fig7.pdf}
\caption{Band structure of the Majorana fermions on a hyperhoneycomb lattice in the zero-flux sector along the high symmetry points in the 3D Brillouin Zone \cite{SETYAWAN2010299}. The special points in the BZ are given by $\Gamma=(0,0,0)$, $Y=\left(0,\frac{2\pi}{b},0\right)$, $T=\left(0,\frac{2\pi}{b},\frac{2\pi}{c}\right)$, $Z=\left(0,0,\frac{2\pi}{c}\right)$, $X=\left(\frac{29\pi}{18a},0,0\right)$ and $A_1=\left(\frac{11\pi}{18a},\frac{2\pi}{b},0\right)$. The fermions with negative energy (bands in black and green) at momentum $\mathbf{k}$ are related to the physical fermions with positive energy (bands in blue and red) at momentum $-\mathbf{k}$ due to particle-hole symmetry.}
\label{fig:ebandstruct}
\end{figure}

\section{Group theoretical analysis  of   lattice vibrations  on the hyperhoneycomb lattice}

\subsection{Counting Raman-Active Phonon Modes

To count the Raman active phonon modes in the hyperhoneycomb lattice, we start by considering the symmetry properties and the number of atoms in the unit cell (12 Oxygens + 4 Iridium + 8 Lithium ions).  The space group for the hyperhoneycomb lattice is the nonsymmorphic $Fddd$ group no. 70. The factor group of the space group over the translation subgroup $Fddd ∖T$ is isomorphic to the point group $D_2h$ \cite{Melnikov2022}. Since the Raman tensors only depend on the point group, we classify the Raman active phonon modes according to the irreps of $D_2h$.  The presence of inversion symmetry divides the phonons into 36 odd infrared active modes ($A_u,B_1u,B_2u,B_3u$) and 36 even Raman active modes ($A_g,B_1g,B_2g,B_3g$).

Next we recall that $D_2h$ can be decomposed as a direct product of $D_2⊗C_i$ since each of the group operations of both groups commute with each other.  For a 3D representation, the generator matrices for $D_2h=D_2⊗C_i$ are given by $𝐱$=diag\{1,-1,-1\}, $𝐲$=diag\{-1,1,-1\} and $𝐢$=diag\{-1,-1,-1\}. These generators can be identified with the group elements $C_2a$, $C_2b $ and inversion respectively.

We choose the origin to be placed at the inversion center (origin choice 2 \cite{ITA}).
With respect to this origin,  we perform the group generator operations $C_2a$, $C_2b $, and inversion and store the transformation of various sites under each symmetry operation in a matrix, which we refer to as the site transformation matrix. This is a $24×24$ matrix (12 Oxygen + 4 Iridium + 8 Lithium ions). The direct product of the 3D generator matrices with the site transformation matrix gives the $72×72$-dimensional matrix representations of the three generators $C_2a$, $C_2b $, and $i$.
Using these matrices for the generators and the group multiplication table, we obtain $72×72$ matrix representations for all 8 symmetry operations
of  $D_2h$ ($E$, $C_2a$, $C_2b$, $C_2c$, $i$,  $σ_bc$, $ σ_ac$, and  $σ_ab$)
in the orthorhombic coordinates $(a,b,c)$. We calculate the character $χ$ of each group element in this naive $72×72$ representation and get :
\[
\begin{array}{c|cccccccc}
& E & C_{2a} & C_{2b} & C_{2c} & i & \sigma_{bc} & \sigma_{ac} & \sigma_{ab}  \\ \hline
\chi & 72 & -12 & 0 & -4 & 0 & 0 & 0 & 0\\
\end{array}
\]
Then, we use the reduction formula to decompose the naive representation into the blocks of irreducible representations:
\begin{equation}\label{eqn:redfor}
\end{equation}
Here $χ^(j)(R)$ is the character of the group element $R$ in the $j^th$ irrep obtained from the character table and $χ(R)$ is the character of group element $R$ obtained in the naive (reducible) representation.
%The reduction formula \eqref{eqn:redfor} enumerates the number of blocks in the $j^{th}$ irrep.
%and  enumerates the number of blocks in $j^{th}$ irrep. %\textcolor{teal}{Thus sentence looks incomplete.}
The irreducible decomposition is then given by
\begin{widetext}
\begin{equation}\label{modect}
\Gamma_\text{Li,Ir,O}=7A_g\oplus8B_{1g}\oplus11B_{2g}\oplus10B_{3g}\oplus 7A_u\oplus8B_{1u}\oplus11B_{2u}\oplus10B_{3u}.
\end{equation}
\end{widetext}
The Raman active modes  of the unit cell are  then
\begin{equation}
\Gamma^{\text R}_\text{Li,Ir,O}=7A_g\oplus8B_{1g}\oplus11B_{2g}\oplus10B_{3g}.
\end{equation}

We note that this factor group analysis of the phonon modes is in agreement with the Raman active modes assigned by Wyckoff positions (WP). The WPs of all ions in $β-Li_2IrO_3$ are \cite{biffin2014unconventional}
\[
\begin{array}{c|ccccc}
& \text{Ir} & \text{Li}_1 & \text{Li}_2 & \text{O}_1 & \text{O}_2 \\ \hline
& 16g & 16g & 16g & 16e & 32h\\
\end{array}
\]
From the selection rule table \cite{bilbaoSAM}, the number of Raman active modes in each irrep channel are
\[\begin{array}{|c||c|c|c|c|}
\hline
\text{WP} & A_g & B_{1g} & B_{2g} & B_{3g} \\
\hline
\hline
16e & 1 & 2 & 2 & 1 \\
\hline
16g & 1 & 1 & 2 & 2 \\
\hline
32h & 3 & 3 & 3 & 3 \\
\hline
\end{array}\]
For example, from the Wyckoff positions, we expect 2 $A_g$ modes corresponding to Li, 1 $A_g$ mode for Ir and 4 $A_g$ modes for O which adds up to the 7 $A_g$ modes. Similarly, for other irreps, we have
\begin{subequations}
\begin{align}
7A_g = 1 (\text{Ir}) + 2 (\text{Li}) + 1 (\text{O}_1) + 3 (\text{O}_2), \nonumber\\
8B_{1g} = 1 (\text{Ir}) + 2 (\text{Li}) + 2 (\text{O}_1) + 3 (\text{O}_2), \nonumber\\
11B_{2g} = 2 (\text{Ir}) + 4 (\text{Li}) + 2 (\text{O}_1) + 3 (\text{O}_2), \nonumber\\
10B_{3g} = 2 (\text{Ir}) + 4 (\text{Li}) + 1 (\text{O}_1) + 3 (\text{O}_2). \nonumber
\end{align}
\end{subequations}
This matches the factor group analysis in Eq.\eqref{modect}.

However, as  explained in Sec. \ref{sec:sp-phCoupling}, in our analysis we only focus on  iridium vibrations. Thus, we restrict the generators to $12×12$-dimensional representation. Using these matrices for the generators and the group multiplication table, we obtain $12×12$ matrix representations for all 8 symmetry operations in the orthorhombic coordinates $(a,b,c)$. We calculate the character $χ$ of each group element in this naive $12×12$ representation-
\[
\begin{array}{c|cccccccc}
& E & C_{2a} & C_{2b} & C_{2c} & i & \sigma_{bc} & \sigma_{ac} & \sigma_{ab}  \\ \hline
\chi & 12 & 0 & 0 & -4 & 0 & 0 & 0 & 0\\
\end{array}
\]
This representation is reduced as
\begin{equation}
\Gamma_\text{Ir}=A_g\oplus B_{1g}\oplus2B_{2g}\oplus2B_{3g}\oplus A_u\oplus B_{1u}\oplus2B_{2u}\oplus2B_{3u}
\end{equation}
and the  Raman active  phonon modes are then classified as
\begin{equation}\label{eqn:IrRaman}
\Gamma^{\text R}_\text{Ir}=A_g\oplus B_{1g}\oplus2B_{2g}\oplus2B_{3g}.
\end{equation}

\subsection{Phonon eigenmode evaluation}\label{subsec:eigenmode}

The aim of this section  is to obtain  the $A_g$ phonon eigenmode, considering vibrations of Ir ions only.   We will do this
by performing a canonical decomposition \cite{Serre1977} of the   representation $Γ$ obtained in Sec.~\ref{subsec:factorGroup}. Specifically, we calculate the coefficients  $u_Γ,i$ in the expansion $u_Γ = ∑_i=1^12 u_Γ, i u_i$, where only vibration of the 4 Iridium atoms is considered.
This is carried out using the projection operator defined as:
\begin{equation}
\mathcal{P}_j=\frac{n_j}{h}\sum_R \, \chi_j(R)^*\,\Gamma_R,
\end{equation}
where $𝒫_j=𝒫_A_g$ is the projection operator for the $A_g$ irrep, $h=8 $ is order of the $D_2h$ group, $n_j=1$ is the dimensionality of the $A_g$ irrep and $χ_j(R)$ are the characters from the character table and $Γ_R$ is the reducible representation for each group operation $R$.  The rank  of $𝒫_A_g$ is 1 since there is  only 1 $A_g$ mode associated with the Ir atoms Eq.~\eqref{eqn:IrRaman}.

The basis vectors of the projection matrix constitute the symmetry-respecting phonon eigenvectors.   These basis vectors are obtained by factorizing the projection operator $𝒫_A_g$ by performing a singular value decomposition (SVD) as $𝒫_A_g=U ·Λ·V^T$. Here $U$ and $V$ are real orthogonal matrices, while $Λ$ is a diagonal matrix with values 0 or 1. Thus, we obtain the basis vectors $u_A_g$ from the columns of the matrix $X_A_g:=U·Λ$. Namely, we have
\begin{equation}
u_{A_g} = \oplus_{i=1}^{12} u_{A_g , i}\, u_i,
\end{equation}
such that $u_A_g=(u^Ir_1_A_g,x,u^Ir_1_A_g,y,u^Ir_1_A_g,zu^Ir_2_A_g,x,u^Ir_2_A_g,y,u^Ir_2_A_g,z,$
The appropriate linear combination of basis vector motivated by DFT results \cite{glamazda2016expt} in orthorhombic coordinates as visualized in Fig. \ref{fig:AgPhonon} is given as: %(and ommiting $\mathbf{q},t$ dependence):
\begin{equation}
\end{equation}

\iffalse
\begin{figure}
\centering
\includegraphics[width=0.95\linewidth]{images/2dos.pdf}
\caption{{\cbl (Don't call it 2DOS - Kinematic Constraints/Allowed processes)} Two-particle density of states (2DOS) as a function of frequency. The particle-hole 2DOS (shown in red) is defined as $\sum_{\bf k}\delta(\Omega-E_{\bf k}^b+E_{\bf k}^a)$. The particle-particle 2DOS (shown in blue) is defined as $\sum_{{\bf k},m,n}\delta(\Omega-E_{\bf k}^m-E_{\bf k}^n)$ where $m,n$ index the two positive bands $a,b$.}
\label{fig:2dos}
\end{figure}
\fi

\section{ Details of the phonon polarization bubble calculation}\label{Appendix:polbub}

In this appendix, we outline the steps  of the derivation of
Eq.~\eqref{eqn:polbubres} for the polarization bubble,  where the total contribution is separated into those from pp- and ph-processes.
The polarization bubble  $ Π(τ)$ from Eq.~(\ref{eqn:polbubexpr}) can be  written as:
\begin{widetext}
\begin{equation}\label{pipptime}
\begin{aligned}
\Pi^{\text{pp}}(\tau)=\frac{\lambda_{A_g}^2}{N}\sum_\mathbf{k}&\bigg[g^a(-\mathbf{k},\tau)g^a(\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{23}\Tilde{\Lambda}_\mathbf{k}^{32} + g^b(-\mathbf{k},\tau)g^b(\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{14}\Tilde{\Lambda}_\mathbf{k}^{41}\\
& + g^a(-\mathbf{k},\tau)g^b(\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{13}\Tilde{\Lambda}_\mathbf{k}^{31} + g^b(-\mathbf{k},\tau)g^a(\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{24}\Tilde{\Lambda}_\mathbf{k}^{42}
\bigg],
\end{aligned}
\end{equation}
\begin{equation}\label{piphtime}
\begin{aligned}
\Pi^{\text{ph}}(\tau)=\frac{\lambda_{A_g}^2}{N}\sum_\mathbf{k}&\bigg[g^a(-\mathbf{k},\tau)\Bar{g}^a(-\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{22}\Tilde{\Lambda}_\mathbf{k}^{22} + g^b(-\mathbf{k},\tau)\Bar{g}^b(-\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{11}\Tilde{\Lambda}_\mathbf{k}^{11}\\
& + g^a(-\mathbf{k},\tau)\Bar{g}^b(-\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{21}\Tilde{\Lambda}_\mathbf{k}^{12} + g^b(-\mathbf{k},\tau)\Bar{g}^a(-\mathbf{k},\tau) \Tilde{\Lambda}_\mathbf{k}^{12}\Tilde{\Lambda}_\mathbf{k}^{21}
\bigg],
\end{aligned}
\end{equation}
\end{widetext}
where the fermionic Green's functions are defined as $g^γ(𝐤,τ)=-⟨T_τγ(τ)γ^†(0)|$⟩ and $\Bar{g}^\gamma(\mathbf{k},\tau)=-\braket{T_\tau \gamma^\dagger(\tau)\gamma(0)}$ with $γ=a, b$ representing the Bogoliubov complex fermion flavor index.
%where $g^\gamma(\mathbf{k},\tau)$, with $\gamma=a,\,b$ representing the complex fermion flavor index, are fermionic  Green's functions:
\iffalse
\begin{equation}
\begin{gathered}
g^\gamma(\mathbf{k},\tau)=-\braket{T_\tau \gamma(\tau)\gamma^\dagger(0)},\\
\Bar{g}^\gamma(\mathbf{k},\tau)=-\braket{T_\tau \gamma^\dagger(\tau)\gamma(0)}.
\end{gathered}
\end{equation}
\fi
For finite temperatures, we perform a Fourier transform to the Matsubara space:
\begin{equation}
\begin{aligned}
g^\gamma(\mathbf{k},i\omega_n)&=\int_0^\beta d\tau\, e^{i\omega_n\tau}\, g^\gamma(\mathbf{k},\tau),\\
g^\gamma(\mathbf{k},\tau)&=\frac{1}{\beta}\sum_{\omega_n} e^{-i\omega_n\tau}\, g^\gamma(\mathbf{k},i\omega_n),
\end{aligned}
\end{equation}
where $g^γ(𝐤,iω_n)= (iω_n-ϵ^γ_𝐤)^-1$ and $g̅^γ(𝐤,iω_n)= (iω_n+ϵ^γ_𝐤)^-1$.

The Fourier transform to  Matsubara space for the polarization bubble is given by:
\begin{equation}
\Pi(i\Omega_m)=\int d\tau \, e^{i\Omega_m\tau}\,\Pi(\tau).
\end{equation}

This Fourier transform results in a convolution of the product of Green's functions appearing in Eq.~\eqref{pipptime} and \eqref{piphtime}. We can evaluate the Matsubara sum for each pair of Green's function as:
\begin{subequations}\label{eqn:matsum}
\allowdisplaybreaks
\begin{align}
\\
\end{align}
\end{subequations}

To obtain real frequencies, we perform an analytic continuation by substituting $iΩ_m→Ω+iδ$.

Finally, we comment on the temperature evolution of pp- and ph-channels of the polarization bubble as shown in Fig. \ref{fig:pitmpplt}. We observe that the pp-channel is dominant compared to ph-channel. This is because an incoming phonon with sufficient energy can decay into a pair of fermions. The kinematic constraints for this process can be satisfied in more ways, especially due to the presence of two physical bands leading to a dominant pp-channel. On the other hand, in a ph-process, a phonon scatters off a fermion and excites it to a higher energy state by transferring the energy of the optical phonon. This process requires a finite population of the fermionic states which only occurs at finite temperatures. Thus, both Re$Π$ and Im$Π$ vanish at zero temperature in the ph-channel.

Furthermore, the renormalized phonon peak  \eqref{eq:broadened_phonon} parameters such as peak position and peak width are controlled by Re$Π$ and Im$Π$ respectively. As noted above, the total contribution obtained by summing over pp- and ph-channels is expected to be pp-dominant. We find that the renormalization to the peak position decreases with temperature as seen from a decrease in Re$Π$ shown by red dashed curve in Fig. \ref{fig:pitmpplt}. This is expected because the spins are not fractionalized at higher temperatures and there are no Majorana fermions for the phonons to couple. The phonon peak width renormalization exhibits a similar decreasing behavior as seen by the decay of Im$Π$ with temperature denoted by red curve in Fig. \ref{fig:pitmpplt}. Thus, with increasing temperature, the effects of the spin-phonon coupling die out.

%\section{Raman intensity calculation}
\section{Magnetic Raman response}\label{app:mag-raman}

In momentum space within the Bogoliubov basis, the magnetic Raman operator (\ref{RamanLF}) can be expressed as follows:
\begin{equation}
\mathcal{R}_{\text{em-s}}=\nu\sum_\mathbf{k}A_\mathbf{k}^\dagger\Tilde{R}_\mathbf{k}A_\mathbf{k},
\end{equation}
where $ν$ is the photon-spin coupling constant,
$R_𝐤=U_𝐤^†R_𝐤 U_𝐤$, and  $U_𝐤$
is the unitary matrix that diagonalizes the spin Hamiltonian  \eqref{eqn:c2a}.
The matrix $R_𝐤$ is given by
\begin{equation}
0 & 0 & \mathcal{B}_{\mathbf{k}} & \mathcal{A}_{\mathbf{k}} \\
0 & 0 & \mathcal{A}_{\mathbf{k}} & \mathcal{C}_{\mathbf{k}} \\
\mathcal{B}^*_{\mathbf{k}} & \mathcal{A}^*_{\mathbf{k}} & 0 & 0 \\
\mathcal{A}^*_{\mathbf{k}} & \mathcal{C}^*_{\mathbf{k}} & 0 & 0
\end{pmatrix},
\end{equation}
with $𝒜_𝐤=-iJ_KP_z$, $ℬ_𝐤=-i J_K(P_x' e^i𝐤·𝐚_2+P_y e^i𝐤·𝐚_1)$, $𝒞_𝐤=-i J_K( P_x+ P_y' e^-i 𝐤·𝐚_3)$,
where we define the polarization factors $P_t=(𝐌_t·ε_in)(𝐌_t·ε_out)$ for the nearest neighbor bonds indexed by $t$.
The Raman intensity is then computed as
\begin{equation}
I_{\text{em-s}}(i\Omega_m)=  -\int_0^\beta  d\tau e^{i\Omega_m \tau}\braket{T_{\tau}\mathcal{R}_{\text{em-s}}(\tau)\mathcal{R}_{\text{em-s}}(0)}.
\end{equation}

\begin{figure}[t]
\centering
\includegraphics[width=0.85\linewidth]{fig8.pdf}
\caption{Temperature evolution of the real and imaginary parts of the polarization bubble, Re$\Pi(\Omega, T)$ and Im$\Pi(\Omega, T)$, at a fixed frequency of $\Omega_{A_g}=24$ meV. The pp- and ph-channels are shown separately.}
\label{fig:pitmpplt}
\end{figure}

We use the same method for calculating this correlator as we did for the polarization bubble in Appendix
\ref{Appendix:polbub}. By expressing the Raman operator in  the Bogoliubov basis, we combine it with the Green's functions describing Majorana fermion's propagators. After performing the Matsubara summation, the result can be written as:

\begin{align}  \label{eqn:Iems-pp}
&I_{\text{em-s}}^{\text{pp}}(i\Omega_m)= \frac{\nu^2}{N}\sum_\mathbf{k} \bigg[P^{gg}_{aa}(i\Omega_m)\Tilde{R}_\mathbf{k}^{23}\Tilde{R}_\mathbf{k}^{32}   +  \\ \nonumber
&P^{gg}_{bb}(i\Omega_m)\Tilde{R}_\mathbf{k}^{14}\Tilde{R}_\mathbf{k}^{41}         +P^{gg}_{ba}(i\Omega_m)\Tilde{R}_\mathbf{k}^{13}\Tilde{R}_\mathbf{k}^{31}  +P^{gg}_{ab}(i\Omega_m)\Tilde{R}_\mathbf{k}^{42}\Tilde{R}_\mathbf{k}^{24}
\bigg],
\end{align}

\begin{align}
\label{eqn:Iems-ph}
&I_{\text{em-s}}^{\text{ph}}(i\Omega_m)= \frac{\nu^2}{N}\sum_\mathbf{k} \Tilde{R}_\mathbf{k}^{43}\Tilde{R}_\mathbf{k}^{34}\bigg[P^{g\bar{g}}_{ab}(i\Omega_m)
+ P^{g\bar{g}}_{ba}(i\Omega_m)
\bigg].
\end{align}
As noted earlier, the contribution of the ph-processes
to the Raman intensity arises only due to the presence of two distinct physical bands, $a$ and $b$ in the Majorana  fermion spectrum  in Kitaev hyperhoneycomb model. This is in stark contrast to the 2D case, where the ph-processes  contribution to the Raman response is not allowed.

Finally, the Raman susceptibility
that characterizes the response of a material to incident light in a Raman scattering experiment at a given frequency
is evaluated as
$χ(Ω)=I_em-s(iΩ_m)|_iΩ_m→Ω+iδ$. The Raman intensity is related to the imaginary part of the Raman susceptibility through the fluctuation-dissipation theorem. The relationship is given by:
\cite{Coleman_2015}:
\begin{equation}
I(\Omega)=-\frac{2}{1-e^{-\beta\Omega}}\;\mathrm{ Im}\chi(\Omega),
\end{equation}
where $β=1/T$
(with $k_B=1$ being the Boltzmann constant). This formula ensures that the intensity is correctly scaled  with temperature according to the Bose-Einstein distribution.\\

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{fig9.pdf}
\caption{Feynman diagrams contributing to the Fano lineshape: (a) and (c) show the contributions from the  ph-channel and (b) and (d) show the contributions from the  pp-channel.  The indices $i,j,l,m \in (1,4)$. The indices $\gamma,\gamma'$ can take values $a$ or $b$, representing the flavor of the complex fermions. The specific values of $\gamma$ and $\gamma'$ are determined by the choice of $i,j,l,m$.}
\label{fig:fanoppph}
\end{figure*}

\section{ Evaluation of the Fano Intensity}\label{App:Fano}
The lowest (first) order diagrams that contribute to the Fano lineshape are given by the following correlation functions:
\begin{subequations}
\begin{align}
I_L^{(1)}(\tau)=\int_0^\beta\, d\tau_1 \,\braket{T_\tau\mathcal{R}^{A_g}_{\text{em-ph}}(\tau)\mathcal{R}^{A_g}_{\text{em-s}}(0)\mathcal{H}^{A_g}_{\text{s-ph}}(\tau_1)}, \\
I_R^{(1)}(\tau)=\int_0^\beta\, d\tau_1 \,\braket{T_\tau\mathcal{R}^{A_g}_{\text{em-s}}(\tau)\mathcal{R}^{A_g}_{\text{em-ph}}(0)\mathcal{H}^{A_g}_{\text{s-ph}}(\tau_1)}.
\end{align}
\end{subequations}
Plugging in the Raman vertices, coupling vertex, and performing Wick's contraction to form propagators from the product of operators yields Eq. \eqref{eqn:IlIr},
\begin{widetext}
\begin{subequations}\label{eqn:IlIr}
\begin{align}
I_L^{(1)}(\tau)&=\int_0^\beta\, d\tau_1 \,\braket{T_\tau\bigg(\mu_{A_g}\mathcal{R}_{A_g}u_{A_g} \bigg)(\tau)\bigg(\nu\sum_\mathbf{k}A_\mathbf{k}^\dagger\Tilde{R}_\mathbf{k}A_\mathbf{k}\bigg)(0)\bigg(\lambda_{A_g}u_{A_g} \sum_\mathbf{k'}A_\mathbf{k'}^\dagger\Tilde{\Lambda}_\mathbf{k'}A_\mathbf{k'}\bigg)(\tau_1)}\\
&=\nu\lambda_{A_g}\mu_{A_g} R^{cc}_{A_g}\int_0^\beta\, d\tau_1 \sum_\mathbf{k}
\braket{T_\tau u_{A_g}(\tau)u_{A_g}(\tau_1)}\braket{T_\tau (A_\mathbf{k}^\dagger)^i(0)A^m_\mathbf{k}(\tau_1)}
\braket{T_\tau A^j_\mathbf{k}(0)(A_\mathbf{k}^\dagger)^l (\tau_1)} \Tilde{R}^{ij}_\mathbf{k}\Tilde{\Lambda}^{lm}_\mathbf{k},\nonumber\\
I_R^{(1)}(\tau)&=\int_0^\beta\, d\tau_1 \,\braket{T_\tau\bigg(\nu\sum_\mathbf{k}A_\mathbf{k}^\dagger\Tilde{R}_\mathbf{k}A_\mathbf{k}\bigg)(\tau)\bigg(\mu_{A_g}\mathcal{R}_{A_g}u_{A_g} \bigg)(0)\bigg(\lambda_{A_g}u_{A_g} \sum_\mathbf{k'}A_\mathbf{k'}^\dagger\Tilde{\Lambda}_\mathbf{k'}A_\mathbf{k'}\bigg)(\tau_1)}\\
&=\nu\lambda_{A_g}\mu_{A_g} R^{cc}_{A_g}\int_0^\beta\, d\tau_1 \sum_\mathbf{k}
\braket{T_\tau (A_\mathbf{k}^\dagger)^i(\tau)A^m_\mathbf{k}(\tau_1)}
\braket{T_\tau A^j_\mathbf{k}(\tau)(A_\mathbf{k}^\dagger)^l (\tau_1)} \braket{T_\tau u_{A_g}(0)u_{A_g}(\tau_1)}\Tilde{R}^{ij}_\mathbf{k}\Tilde{\Lambda}^{lm}_\mathbf{k}.\nonumber
\end{align}
\end{subequations}
\end{widetext}

By performing a Fourier transform to Matsubara frequency space and subsequently carrying out the Matsubara sums, we derive the final expressions for the Fano intensity, as given in Eq.~\eqref{eqn:fano}. %\textcolor{teal}{incorrect ref?} in the main text.
The left and right Fano diagrams, with the ph- and pp-channels explicitly shown, are presented in Fig.~\ref{fig:fanoppph}. Their explicit expressions are
\begin{widetext}
\begin{subequations}\label{eqn:fano}
\begin{equation}\label{eqn:ilpp}
\begin{aligned}
I_L^{\text{pp}}(i\Omega_m)=-\frac{\lambda_{A_g}\nu\mu_{A_g}}{N}\sum_\mathbf{k} \mathcal{D}(i\Omega_m)\bigg[P_{aa}^{gg}(i\Omega_m)\Tilde{\Lambda}_\mathbf{k}^{23}\Tilde{R}_\mathbf{k}^{32} + P_{bb}^{gg}(i\Omega_m)\Tilde{\Lambda}_\mathbf{k}^{14}\Tilde{R}_\mathbf{k}^{41}
+ P_{ba}^{gg}(i\Omega_m) \Tilde{\Lambda}_\mathbf{k}^{13}\Tilde{R}_\mathbf{k}^{31} + P_{ab}^{gg}(i\Omega_m)\Tilde{\Lambda}_\mathbf{k}^{24}\Tilde{R}_\mathbf{k}^{42}
\bigg],
\end{aligned}
\end{equation}
\begin{equation}\label{eqn:irpp}
\begin{aligned}
I_R^{\text{pp}}(i\Omega_m)=-\frac{\lambda_{A_g}\nu\mu_{A_g}}{N}\sum_\mathbf{k} \mathcal{D}(i\Omega_m)\bigg[P_{aa}^{gg}(i\Omega_m)\Tilde{R}_\mathbf{k}^{23}\Tilde{\Lambda}_\mathbf{k}^{32} + P_{bb}^{gg}(i\Omega_m)\Tilde{R}_\mathbf{k}^{14}\Tilde{\Lambda}_\mathbf{k}^{41}
+ P_{ba}^{gg}(i\Omega_m) \Tilde{R}_\mathbf{k}^{13}\Tilde{\Lambda}_\mathbf{k}^{31} + P_{ab}^{gg}(i\Omega_m)\Tilde{R}_\mathbf{k}^{24}\Tilde{\Lambda}_\mathbf{k}^{42}
\bigg],
\end{aligned}
\end{equation}
\begin{equation}\label{eqn:ilph}
\begin{aligned}
I_L^{\text{ph}}(i\Omega_m)=-\frac{\lambda_{A_g}\nu\mu_{A_g}}{N}\sum_\mathbf{k} \mathcal{D}(i\Omega_m)\bigg[P_{ab}^{\Bar{g}g}(i\Omega_m)\Tilde{\Lambda}_\mathbf{k}^{34}\Tilde{R}_\mathbf{k}^{43} + P_{ba}^{g\bar{g}}(i\Omega_m)\Tilde{R}_\mathbf{k}^{12}\Tilde{\Lambda}_\mathbf{k}^{21}
\bigg],
\end{aligned}
\end{equation}
\begin{equation}\label{eqn:irph}
I_R^{\text{ph}}(i\Omega_m)=-\frac{\lambda_{A_g}\nu\mu_{A_g}}{N}\sum_\mathbf{k} \mathcal{D}(i\Omega_m) P_{ba}^{g\Bar{g}}(i\Omega_m) \bigg[\Tilde{R}_\mathbf{k}^{12}\Tilde{\Lambda}_\mathbf{k}^{21} + \Tilde{R}_\mathbf{k}^{34}\Tilde{\Lambda}_\mathbf{k}^{43}
\bigg],
\end{equation}
\end{subequations}
\end{widetext}
where $𝒟(iΩ_m)$  is the phonon propagator  renormalized via the  usual Dyson equation.

\bibliography{refs.bib}
\bibliographystyle{apsrev4-1}
\end{document}