# Kondo destruction and fixed-point annihilation in a Bose-Fermi Kondo model

Haoyu Hu<EMAIL_ADDRESS>Department of Physics & Astronomy, Rice Center for
Quantum Materials, Rice University, Houston, Texas 77005, USA  Qimiao Si
<EMAIL_ADDRESS>Department of Physics & Astronomy, Rice Center for Quantum
Materials, Rice University, Houston, Texas 77005, USA

###### Abstract

Quantum criticality that goes beyond the Landau framework of order-parameter
fluctuations is playing a central role in elucidating the behavior of strange
metals. A prominent case appears in Kondo lattice systems, which have been
extensively analysed in terms of an effective Bose-Fermi Kondo model. Here, a
spin is simultaneously coupled to conduction electron bands and gapless vector
bosons that represent magnetic fluctuations. The Bose-Fermi Kondo model
features interacting fixed points of Kondo destruction with such properties as
dynamical Planckian ($\hbar\omega/k_{\rm B}T$) scaling and loss of
quasiparticles. Here we carry out a renormalization-group analysis of the
model with spin isotropy and identify pair-wise annihilations of the fixed
points as the spectrum of the bosonic bath evolves. Our analysis not only
provides an essentially complete understanding of the previous numerical
results of an SU(2)-symmetric model, but also reveals a surprising feature of
sequential fixed-point annihilation. Our results lay the foundation for the
understanding of quantum criticality in spin-isotropic heavy-fermion metals as
well as in doped Mott-Hubbard systems.

Introduction:  Quantum criticality appears in a variety of strongly correlated
metallic systems Special issue (2013); Coleman and Schofield (2005); Paschen
and Si (2021); Kirchner  _et al._ (2020); Sachdev (1999), with the
antiferromagnetic (AF) case being prototypical. It provides a setting for
strange metals with a complete loss of quasiparticles and with a dynamical
Planckian ($\hbar\omega/k_{\rm B}T$) scaling Si  _et al._ (2001); Coleman  _et
al._ (2001); Senthil  _et al._ (2004). Such properties are to be contrasted
with the Landau description of metallic quantum criticality. In the Landau
framework, phases of matter are identifieded by order parameters, which
specify the spontaneous breaking of global symmetries, and quantum criticality
is described by the spatial and temporal fluctuations of the order parameter:
A continuous AF to paramagnetic phase transition at $T=0$ is characterized by
the slow fluctuations of the staggered magnetization Hertz (1976); Millis
(1993); at such a spin-density-wave QCP, quasiparticles remain and the
underlying fixed point is Gaussian, leading to a violation of dynamical
Planckian scaling. In AF heavy-fermion metals, the beyond-Landau description
Si  _et al._ (2001); Coleman  _et al._ (2001); Senthil  _et al._ (2004)
features the added critical ingredient in the form of Kondo destruction.
Experiments in many heavy fermion strange metals are in strong support for the
beyond-Landau nature of the quantum critical points Paschen and Si (2021);
Kirchner  _et al._ (2020); Paschen  _et al._ (2004); Gegenwart  _et al._
(2007); Friedemann  _et al._ (2010); Shishido  _et al._ (2005); Park  _et al._
(2006); Knebel  _et al._ (2008); Custers  _et al._ (2012); Martelli  _et al._
(2019); Schröder  _et al._ (2000); Prochaska  _et al._ (2020); Chen  _et al._
(2022).

Given the broad current interest in the strange metal behavior Keimer and
Moore (2017); Paschen and Si (2021); Phillips  _et al._ (2022), it is
opportune to elucidate its physics by expanding our understanding of the Kondo
destruction quantum criticality. The Bose-Fermi Kondo (BFK) model has been
extensively recognized as an important case study in this context, given that
it is an effective description Si  _et al._ (2014) of the Kondo lattice
Hamiltonian in the extended dynamical mean-field theory (EDMFT) Si and Smith
(1996a); Smith and Si (2000); Chitra and Kotliar (2000). The model has been
treated by renormalization group (RG) method within an $\epsilon$-expansion
(see below for the definition of $\epsilon$) Si and Smith (1996b); Smith and
Si (1999); Sengupta (2000); Si  _et al._ (2001, 2003); Zhu and Si (2002);
Zaránd and Demler (2002). Various numerical means ranging from Monte-Carlo to
numerical renormalization group methods Grempel and Si (2003); Zhu  _et al._
(2003); Glossop and Ingersent (2007); Zhu  _et al._ (2007) have been used to
analyse the Ising-anisotropic model. For the SU(2)-symmetric BFK model,
numerical studies have only recently become possible: Ref. Cai and Si, 2019
developed a continuous-time Monte Carlo method for this case and identified an
unexpectedly large number (eight in total) of fixed points for small (but
positive) $\epsilon$, which are pair-wise annihilated as $\epsilon$ is
increased. Understanding these surprising results is important not only for
the beyond-Landau quantum criticality and strange metal physics but also for
the decoherence physics of spin qubits.

In this work, we carry out analytical calculations to identify all the fixed
points of the SU(2)-symmetric BFK model seen in the numerical simulations. Our
calculation utilizes an RG method that is based on a large-$S$ expansion (with
$S$ denoting the spin size) Cuomo  _et al._ (2022). Our results not only
capture the phenomena of pair-wise fixed-point annihilations as a function of
increasing $\epsilon$, as suggested by the numerical results, but also uncover
a hitherto unsuspected result that such annihilations develop in sequence. The
BFK model can be understood as an effective $0+1$ dimensional nonlinear
$\sigma$ model, but with not only a topological term but also an extra Kondo
coupling to gapless fermions. As we show later, the Kondo coupling affects the
topological term. The interplay among the Kondo interaction, topological term,
and bosonic coupling leads to a rich phase diagram with interacting fixed
points that annihilate in sequence. Our results are important for the
understanding of quantum criticality in correlated systems such as the physics
of Kondo destruction in heavy fermion metals and doped Mott-Hubbard systems.

Model:  We focus on the following Hamiltonian for the BFK model:

$\displaystyle H$ $\displaystyle=$ $\displaystyle
J_{K}\sum_{\mu,a}S^{\mu}s^{\mu}_{c,a}(\bm{r_{0}})+g\sum_{\mu}S^{\mu}\phi^{\mu}(\bm{r_{0}})$
(1) $\displaystyle+$
$\displaystyle\sum_{\bm{k},\sigma,a}\epsilon_{\bm{k}}c_{\bm{k},\sigma,a}^{\dagger}c_{\bm{k},\sigma,a}+\frac{1}{2}\sum_{\bm{p},\mu}E_{\bm{p}}\phi_{-\bm{p}}^{\mu}\phi_{\bm{p}}^{\mu}\,.$

Here, $S^{\mu}(\mu\in\\{x,y,z\\})$ denotes a spin-$S$ impurity at site
$\bm{r_{0}}=\bm{0}$. $c_{\bm{k},\sigma,a}^{\dagger}$ creates a conduction
electron with momentum $\bm{k}$, spin $\sigma$, orbital $a$ (where
$a\in\\{1,2,..,K\\}$) and dispersion $\epsilon_{\bm{k}}$. We consider a
generic filling of the conduction electrons with a Fermi surface in the non-
interacting limit, which gives a constant conduction electron density of state
near the Fermi energy:
$\rho_{c}(\epsilon)=\sum_{\bm{k}}\delta(\epsilon-\epsilon_{\bm{k}})\sim\rho_{0}$.
The spin operator of the conduction electron with orbital $a$ at site
$\bm{r_{0}}$ is defined as
$s_{c,a}^{\mu}(\bm{r_{0}})=c_{\bm{r_{0}},a}^{\dagger}\frac{\sigma^{\mu}}{2}c_{\bm{r_{0}},a}$,
where $\sigma^{x,y,z}$ are the Pauli matrices. It couples to the impurity at
the same site via Kondo coupling $J_{K}$. $\phi_{\bm{p}}^{\mu}$ describes a
free $O(3)$ bosonic field with momentum $\bm{p}$ and kinetic energy
$E_{\bm{p}}=\bm{p}^{2}$. The spectral functions of $\phi$ fields are taken to
be $\rho_{b}(\omega)=\sum_{\bm{p}}\delta(\omega-
E_{\bm{p}})\propto|\omega|^{1-\epsilon}$. $g$ is the coupling strength between
bosonic field at site $\bm{r}_{0}$ and impurity spin. To perform a large $S$
expansion, we also rescale the fields and coupling constants as
$\phi^{\mu}\rightarrow\sqrt{S}\phi^{\mu}$, $g\rightarrow S^{-1/2}g$,
$J_{K}\rightarrow 2S^{-1}J_{K}$. All the terms of the Hamiltonian are then of
order $S$. Throughout the work, we focus on the case of perfect screening with
$\kappa=K/S=2$.

RG analysis and beta functions: We now carry out an RG analysis of the model
defined in Eq. 1. To do so, we first calculate the propagators of the fermions
and bosons using the large $S$ expansion.

We note that these propagators can be expressed as functionals of the
$n$-point correlation functions of $S^{\mu}$ [see Supplementary Material (SM)
A SM ]. Thus, it’s sufficient to consider the effective theory of $S^{\mu}$
and its correlation functions. The effective action of $S^{\mu}$ reduces to

$\displaystyle S_{eff}$ $\displaystyle=$ $\displaystyle-
iSS_{B}-\frac{g^{2}}{2}\int_{\tau,\tau^{\prime}}\sum_{\tau,\tau^{\prime}}S^{\mu}(\tau)G_{\phi}(\tau-\tau^{\prime})S^{\mu}(\tau^{\prime})$
$\displaystyle-K\text{Tr}\log\bigg{[}\delta_{\tau,\tau^{\prime}}\bigg{(}\delta_{k,k^{\prime}}(\partial_{\tau^{\prime}}+\epsilon_{\bm{k}})+J_{K}\sum_{\mu}S^{\mu}(\tau)\frac{\sigma^{\mu}}{2}\bigg{)}\bigg{]}\,,$

where $S_{B}$ describes the Berry phase for the spin, the kernel
$G_{\phi}(\tau)\propto|\tau-\tau^{\prime}|^{-(2-\epsilon)}$ is obtained by
integrating out the bosonic fields and the last line represents the
contribution from integrating out the conduction electrons. We then introduce
the spinor representation of the impurity spin
$S^{\mu}=Sz_{a}^{\dagger}\frac{\sigma_{ab}^{\mu}}{2}z_{b}$ with a constraint,
$|z_{1}|^{2}+|z_{2}|^{2}=2$ Auerbach (2012). We expand the action around the
saddle point solution $z=z_{0}=\text{const}$ and integrate out zero-modes to
enforce the $SU(2)$ symmetry Zinn-Justin (2007). Without loss of generality,
we consider $z=z_{0}+\delta z$ with

$\displaystyle
z_{0}=\sqrt{2}\begin{bmatrix}1&0\end{bmatrix}\,,~{}~{}~{}~{}~{}~{}~{}~{}\delta
z=\begin{bmatrix}-\frac{|\chi|^{2}}{\sqrt{2}S}&\frac{i\chi^{\dagger}}{\sqrt{S}}\end{bmatrix}\,.$
(3)

Here, the $\chi$ fields parametrize the $1/S^{0}$ fluctuations around the
saddle-point solution $z_{0}$ Auerbach (2012). The effective action at order
$S^{0}$ is

$\displaystyle S_{\chi}=\int
G^{-1}_{\chi}(i\Omega)|\chi(i\Omega)|^{2}\frac{d\Omega}{2\pi}\,,$ (4)

with the propagator (see SM A SM ):

$\displaystyle G_{\chi}^{-1}(i\Omega)$ $\displaystyle=$ $\displaystyle
i\Omega\bigg{[}1+p(J|\Omega|^{-r})\bigg{]}+|\Omega|\bigg{[}q(J|\Omega|^{-r})+h|\Omega|^{-\epsilon}\bigg{]}\,.$

Here, the functions $q(x)$ and $p(x)$ are defined as

$\displaystyle q(x)=2\kappa x^{2}\,/\left[\pi\,(1+x^{2})^{2}\right]\,,$
$\displaystyle
p(x)=(\kappa/\pi)[x-x^{3}-(1+x^{2})^{2}\arctan(x)]/(1+x^{2})^{2}\,,$

where $p(x)\sim-8\kappa x^{3}/(3\pi)$ at small $x$, and
$p(x\rightarrow\infty)\rightarrow-\kappa/2$. In the derivation, we regularize
the conduction electron density of state to be
$\rho_{c}(\epsilon)=\rho_{0}|\epsilon|^{-r}$. The exponent $r$ is introduced
to perform a minimal-subtraction RG study and will be set to zero at the final
step of the calculation Zhu and Si (2002). Other constants in Eq.
LABEL:eq:prop are $h\propto g^{2}$, $J=\pi J_{K}\rho_{0}$.

With the effective theory of $\chi$, which describes the $S^{0}$ fluctuations,
we’re able to calculate the correlation functions of the fermions and bosons
(see SM A SM ). We then require that the poles in these correlation functions
are minimally removed Zhu and Si (2002); Zinn-Justin (2007), which gives the
following beta functions (see SM A SM ):

$\displaystyle\beta_{J}$ $\displaystyle=$
$\displaystyle\frac{J-J^{3}}{2S\pi(1+J^{2})}F(J)+\frac{J^{2}}{S\pi(1+J^{2})}H(J)$
$\displaystyle+\bigg{[}\frac{-(1-J^{2})Jh+2iJ^{2}h}{2S(1+J^{2})\pi
c(J)(c(J)+h)}+\text{h.c.}\bigg{]}+o(\frac{1}{S})\,,$ $\displaystyle\beta_{h}$
$\displaystyle=$
$\displaystyle-\frac{1}{S}\tilde{\epsilon}h+\frac{h}{S\pi}F(J)-\frac{1}{S\pi}\frac{(c(J)+c(J)^{*})h}{|c(J)|^{2}}$
(6)
$\displaystyle-\frac{1}{S\pi}\frac{-2h^{2}-(c(J)+c(J)^{*})h}{(c(J)+h)(c(J)^{*}+h)}+o(\frac{1}{S})\,.$

Here,

$\displaystyle F(j)$ $\displaystyle=$
$\displaystyle{2q(j)}/{\left[(1+p(j))^{2}+q(j)^{2}\right]}$ $\displaystyle
H(j)$ $\displaystyle=$
$\displaystyle{-2(1+p(j))}/{\left[(1+p(j))^{2}+q(j)^{2}\right]}$
$\displaystyle c(J)$ $\displaystyle=$ $\displaystyle i(1+p(J))+q(J),$ (7)

and we rescale $\epsilon$ as $\epsilon=\tilde{\epsilon}/S$. All terms in the
beta functions at the order of $1/S$ have been captured. By treating $1/S$ as
the only small parameter in the beta function, we’re able to determine the
phase diagram even for large $J$ and $h$.

Figure 1: Schematic renormalization group flows for
$\tilde{\epsilon}<\tilde{\epsilon}_{1}$ (a),
$\tilde{\epsilon}_{1}<\tilde{\epsilon}<\tilde{\epsilon}_{2}$ (b) and
$\tilde{\epsilon}_{2}<\tilde{\epsilon}$ (c). The squares label stabled fixed
points, the circles mark unstable fixed points and the wide line denotes a
line of stable fixed points at $h\rightarrow\infty$. See SM C for the RG flow
from an explicit numerical evaluations of the beta functions.

RG flow:  We are now in a position to analyze the RG flow. Note that there is
always a trivial unstable fixed point at $J=h=0$, which we will not discuss
further. In addition, at $J=0$, the model becomes a Bose-Kondo model and our
results are consistent with previous works Cuomo  _et al._ (2022); Nahum
(2022). After numerical evaluations of Eq. 6 at perfect screening $\kappa=2$,
we identify three different types of phase diagrams corresponding to three
regions of $\tilde{\epsilon}$: $\tilde{\epsilon}<\tilde{\epsilon}_{1}$,
$\tilde{\epsilon}_{1}<\tilde{\epsilon}<\tilde{\epsilon}_{2}$,
$\tilde{\epsilon}_{2}<\tilde{\epsilon}$, with
$\tilde{\epsilon}_{1}=1/\pi,\tilde{\epsilon}_{2}\simeq 0.36$. The schematic RG
flows are shown in Fig. 1. There are in total seven types of non-trivial fixed
points with their properties – including their labels – listed in Tab. 1.

We first analyze the RG flow at $\tilde{\epsilon}<\tilde{\epsilon}_{1}$ as
shown in Fig. 1 (a). At $h=0$, the model is equivalent to a standard Kondo
impurity model and the behavior of $\beta_{J}$ at small $J$ (see SM C SM ) is
consistent from what has been found via perturbation theory Affleck (1995).
Numerically, we find $J$ is always relevant and the system flows to a strong
coupling fixed point at $J=\infty$. We call it Kondo or K fixed point; the
bosonic coupling $h$ is irrelevant at this fixed point. For $J=0$, there are
two stable fixed points at
$h=(1-\sqrt{1-\tilde{\epsilon}^{2}})/{\tilde{\epsilon}},\infty$ and one
critical point at $h=(1+\sqrt{1-\tilde{\epsilon}^{2}})/{\tilde{\epsilon}}$.
Due to their local moment nature, we label the two stable fixed points as L
fixed point (at small $h$) and L′ fixed point (at infinity $h$). We call the
critical point that controls the quantum phase transition between L and L′ as
LC. After introducing the Kondo coupling, we find $J$ is irrelevant at both L
and LC, but marginal at L′. The marginal behavior extends to finite $J$ which
yields a line of fixed points at $h=\infty$ (see SM C SM ).

| Position | Stable? | Region | Color
---|---|---|---|---
L | $J=0,h=\frac{1-\sqrt{1-\tilde{\epsilon}^{2}}}{\tilde{\epsilon}}$ | Stable | $\tilde{\epsilon}<\tilde{\epsilon}_{1}$ | Grey
L′ | $h=\infty$ | Stable |  | Yellow
K | $J=\infty,h=0$ | Stable |  | Purple
LC | $J=0,h=\frac{1+\sqrt{1-\tilde{\epsilon}^{2}}}{\tilde{\epsilon}}$ | Unstable | $\tilde{\epsilon}<\tilde{\epsilon}_{1}$ | Green
KD | $J\neq 0,h\neq 0$ | Unstable | $\tilde{\epsilon}<\tilde{\epsilon}_{2}$ | Red
KD′ | $J=\infty,h=\frac{2}{\tilde{\epsilon}\pi}$ | Unstable |  | Blue
C | $J\neq 0,h\neq 0$ | Unstable | $\tilde{\epsilon}<\tilde{\epsilon}_{2}$ | Brown
Table 1: List of fixed points with their positions in the phase diagram,
stability, region of existence and color in Fig. 1

Away from $J=0$, $h=0$ and $h=\infty$ lines, there are three unstable fixed
points. We label them as KD, KD′, and C. Both KD and KD′ describe the Kondo
destruction quantum phase transition with one relevant direction and one
irrelevant direction. A third fixed point C with all the directions being
relevant separates KD and KD′ in the RG flow. KD is located at smaller $J,h$
and controls the Kondo destruction quantum phase transition between the local
momentum phase L and Kondo-screened phase K. The quantum phase transition
between L′ and K is described by the fixed point KD′ at
$(h_{c}=\frac{2}{\tilde{\epsilon}\pi},J=\infty)$ (see SM C SM ).

The topology of RG flow and the existence of fixed points change as we
increase the value of $\tilde{\epsilon}$. Remarkably, we find a two-step
sequence of fixed-point annihilation. The first annihilation occurs between L
and LC′ at $\tilde{\epsilon}=\tilde{\epsilon}_{1}=1/\pi$ and the second fixed-
point annihilation occurs between C and KD at
$\tilde{\epsilon}=\tilde{\epsilon}_{2}\simeq 0.36$. We show the RG flow after
first annihilation in Fig. 1 (b) and the flow after second annihilation in
Fig. 1 (c).

Figure 2: Zeros of $\beta_{J}$ and $\beta_{h}$ at
$\tilde{\epsilon}<\tilde{\epsilon}_{1}$,
$\tilde{\epsilon}=\tilde{\epsilon}_{1}$,
$\tilde{\epsilon}_{1}<\tilde{\epsilon}<\tilde{\epsilon}_{2}$,
$\tilde{\epsilon}=\tilde{\epsilon}_{2}$ and
$\tilde{\epsilon}>\tilde{\epsilon}_{2}$. $\beta_{h}=0$ line have two crossing
points with $J=0$ at $\tilde{\epsilon}<\tilde{\epsilon}_{1}$, which describe L
and LC fixed points respectively. Two crossing points merge into one at
$\tilde{\epsilon}=\tilde{\epsilon}_{1}$ and disappear at
$\tilde{\epsilon}>\tilde{\epsilon}_{1}$. $\beta_{h}=0$ and $\beta_{J}=0$ lines
have two crossing points at $\tilde{\epsilon}<\tilde{\epsilon}_{2}$, which
represent KD and C fixed points respectively. They merge into one at
$\tilde{\epsilon}=\tilde{\epsilon}_{2}$ and disappear at
$\tilde{\epsilon}>\tilde{\epsilon}_{2}$.

Fixed-point annihilation: We next turn to a detailed analysis of the fixed-
point annihilation. The annihilation between two fixed points respectively
with $n$ and $n+1$ relevant directions is generically allowed by the topology
of the RG flow. In the BFK model, both $n=0$ and $n=1$ types of annihilation
happen. To be specific, we consider the zeros of beta functions. The shape of
$\beta_{J}=0$ line and $\beta_{h}=0$ line at various $\tilde{\epsilon}$ is
shown in Fig. 2. At $\tilde{\epsilon}<\tilde{\epsilon}_{1}$, $\beta_{h}=0$
line have two crossings points with $J=0$ axis, which label the positions of L
(at a smaller $h$) and LC (at a larger $h$). With increasing
$\tilde{\epsilon}$, the crossing points get closer and merge into one at
$\tilde{\epsilon}=\tilde{\epsilon}_{1}$, signaling the first fixed-point
annihilation between L and LC.

The two crossing points between $\beta_{J}=0$ and $\beta_{h}=0$ lines
represent KD (at a smaller $h$) and C (at a larger $h$) fixed points
respectively. Their annihilation occurs at
$\tilde{\epsilon}=\tilde{\epsilon}_{2}$, where two curves have a single
touching point. Upon further increasing $\tilde{\epsilon}$, these fixed points
disappear.

Role of Berry phase:  We now discuss the role of Berry phase and its interplay
with Kondo coupling. From Eq. LABEL:eq:prop, we see that the Kondo coupling
contributes an $i\Omega p(J)$ term to $G_{\chi}^{-1}$ in the $r=0$ limit. It
reduces the Berry phase term, $i\Omega$, at finite $J$ and fully removes it at
$J=\infty$ and $K=2S$. The exact cancellation reveals the nature of Kondo
effect which can also be seen by directly solving the problem at $J=\infty$ as
shown in SM D SM . It also suggests that, at $J=\infty$, the system should
always stay at the Kondo-screened phase. Thus, we expect the critical point
KD′, located at $J=\infty,h=2/(\tilde{\epsilon\pi})$, will be moved to a
finite $J$ when the next-order corrections in $1/S$ are included in the RG
analysis.

Figure 3: RG flow based on numerical simulations at $S=1/2$ (adapted from Ref.
Cai and Si (2019)), where $\tilde{\epsilon}^{*}\simeq 0.265$

Comparison with numerical results: The BFK model is equivalent to the
numerically studied Bose-Fermi Anderson model with a power-law bosonic bath
spectral functions $\rho_{b}(\omega)=|\omega|^{1-\epsilon}$ when the Hubbard
interaction $U$ is sufficiently large Cai and Si (2019). At $S=1/2$, two types
of RG flow have been observed in the numerical simulations as shown in Fig. 3
(adapted from Ref. Cai and Si (2019)). The first one, at
$\tilde{\epsilon}<\tilde{\epsilon}^{*}\simeq 0.265$, is comparable with our
analytical results at $\tilde{\epsilon}<\tilde{\epsilon}_{1}$. The other one,
at $\tilde{\epsilon}>\tilde{\epsilon}^{*}$, is consistent with our analytical
RG flow at $\tilde{\epsilon}>\tilde{\epsilon}_{2}$. We note that our
analytical analysis has identified all the fixed points observed numerically
and the values of $\tilde{\epsilon}_{1,2}$ that we have derived turns out to
be comparable to the numerical results. However, in further details there are
several distinctions in our analytical results. First, from the analytical
calculations the KD′ fixed point is located at infinite $J$, while it is at a
finite $J$ numerically. As we explained before, the next-order corrections (in
$1/S$) are expected to move the position of KD′ to finite $J$ and we leave
their analysis for a future study. Second, the analytical calculations find a
line of fixed points at infinite $h$ instead of a single point as suggested by
numerical simulations. In practice, it’s hard to conclude whether L′ phase
corresponds to a single fixed point or a line of fixed points numerically, and
it’s also possible that the next order (in $1/S$) corrections would shrink the
line of fixed points to a single fixed point in our analytical approaches. We
leave the precise reconcilation of this issue for future studies. Finally,
numerical calculations fail to find RG flows corresponding to
$\tilde{\epsilon}_{1}<\tilde{\epsilon}<\tilde{\epsilon}_{2}$. The likely
reason for this discrepancy is that numerical calculations have so far focused
on only a small number of $\tilde{\epsilon}$ and may have missed this narrow
region of $\tilde{\epsilon}$.

Discussion: Several remarks are in order. First, our results have some general
connections with the dimensional hierarchy of nonlinear $\sigma$ model with
topological terms Nahum  _et al._ (2015); Ma and Wang (2020); Abanov and
Wiegmann (2000). An important distinction of our model is the additional
coupling to gapless fermions, which normalize the effect of the topological
term and leads to a richer phase diagram.

Second, our model and results are in general relevant to various strongly
correlated systems, including heavy fermion metals and doped Mott-Hubbard
systems. In these systems, the electrons tend to be localized and form local
moments due to the strong Coulomb interactions that they experience. The local
moment then interacts with the remaining itinerant electrons and also
collective bosonic excitations, the effects of which are captured by the BFK
model. Importantly, the fixed point KD′ has a local spin susceptibility
$1/|\tau|^{\epsilon}$ Cai and Si (2019), which is essential for realizing the
EDMFT solutions for the beyond-Landau quantum criticality Si  _et al._ (2014).
Thus, the understanding we have reached for the Kondo destruction physics of
the SU(2) BFK model is crucially important for understanding the quantum
criticality of magnetic heavy-fermion metals Special issue (2013); Coleman and
Schofield (2005); Paschen and Si (2021); Kirchner  _et al._ (2020).
Generalizations to including additional local degrees of freedom, such as
multipoles Martelli  _et al._ (2019); Liu  _et al._ (2021); Han  _et al._
(2022), may also be considered. Finally, the BFK model also represents an
effective description of the doped Mott systems in the form of $t$-$J$-$U$
Hubbard models Smith and Si (2000). As such, our results may also play an
important role in the understanding of such systems Fang  _et al._ (2022);
Badoux  _et al._ (2016); Chowdhury  _et al._ (2021).

Conclusions: We have carried out a renormalization group analysis of the spin-
isotropic Bose-Fermi Kondo model. Our RG analysis determines a set of
interacting fixed points with dynamical Planckian scaling and loss of
quasiparticles. Our results not only provide an understanding of the numerical
results on the systematic evolution of the fixed points, but also reveals an
unexpected phenomenon of sequential fixed-point annihilations. Our results lay
the foundation for the understanding of quantum criticality in both spin-
isotropic heavy-fermion metals and doped Mott-Hubbard systems.

We thank A. Cai and S. Paschen for useful discussions. This work has primarily
been supported by the NSF Grant No. DMR-2220603 and by the Robert A. Welch
Foundation Grant No. C-1411. The work of Q.S. was performed in part at the
Aspen Center for Physics, which is supported by the NSF grant No. PHY-1607611.

## References

  * Special issue (2013) Special issue, Phys. Status Solidi 250, 417 (2013).
  * Coleman and Schofield (2005) P. Coleman and A. J. Schofield, Nature 433, 226 (2005).
  * Paschen and Si (2021) S. Paschen and Q. Si, Nat. Rev. Phys. 3, 9 (2021).
  * Kirchner  _et al._ (2020) S. Kirchner, S. Paschen, Q. Chen, S. Wirth, D. Feng, J. D. Thompson, and Q. Si, Rev. Mod. Phys. 92, 011002 (2020).
  * Sachdev (1999) S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, Cambridge, 1999).
  * Si  _et al._ (2001) Q. Si, S. Rabello, K. Ingersent, and J. L. Smith, Nature 413, 804 (2001).
  * Coleman  _et al._ (2001) P. Coleman, C. Pépin, Q. Si, and R. Ramazashvili, J. Phys. Cond. Matt. 13, R723 (2001).
  * Senthil  _et al._ (2004) T. Senthil, M. Vojta, and S. Sachdev, Phys. Rev. B 69, 035111 (2004).
  * Hertz (1976) J. Hertz, Phys. Rev. B 14, 1165 (1976).
  * Millis (1993) A. J. Millis, Phys. Rev. B 48, 7183 (1993).
  * Paschen  _et al._ (2004) S. Paschen, T. Lühmann, S. Wirth, P. Gegenwart, O. Trovarelli, C. Geibel, F. Steglich, P. Coleman, and Q. Si, Nature 432, 881 (2004).
  * Gegenwart  _et al._ (2007) P. Gegenwart, T. Westerkamp, C. Krellner, Y. Tokiwa, S. Paschen, C. Geibel, F. Steglich, E. Abrahams, and Q. Si, Science 315, 1049 (2007).
  * Friedemann  _et al._ (2010) S. Friedemann, N. Oeschler, S. Wirth, C. Krellner, C. Geibel, F. Steglich, S. Paschen, S. Kirchner, and Q. Si, Proceedings of the National Academy of Sciences 107, 14547 (2010).
  * Shishido  _et al._ (2005) H. Shishido, R. Settai, H. Harima, and Y. Ōnuki, J. Phys. Soc. Jpn. 74, 1103 (2005).
  * Park  _et al._ (2006) T. Park, F. Ronning, H. Q. Yuan, M. B. Salamon, R. Movshovich, J. L. Sarrao, and J. D. Thompson, Nature 440, 65 (2006).
  * Knebel  _et al._ (2008) G. Knebel, D. Aoki, J.-P. Brison, and J. Flouquet, J. Phys. Soc. Jpn. 77, 114704 (2008).
  * Custers  _et al._ (2012) J. Custers, K. A. Lorenzer, M. Müller, A. Prokofiev, A. Sidorenko, H. Winkler, A. M. Strydom, Y. Shimura, T. Sakakibara, R. Yu, Q. Si, and S. Paschen, Nat. Mater. 11, 189 (2012).
  * Martelli  _et al._ (2019) V. Martelli, A. Cai, E. M. Nica, M. Taupin, A. Prokofiev, C.-C. Liu, H.-H. Lai, R. Yu, K. Ingersent, R. Küchler, A. M. Strydom, D. Geiger, J. Haenel, J. Larrea, Q. Si, and S. Paschen, Proc. Natl. Acad. Sci. U.S.A. 116, 17701 (2019).
  * Schröder  _et al._ (2000) A. Schröder, G. Aeppli, R. Coldea, M. Adams, O. Stockert, H. v. Löhneysen, E. Bucher, R. Ramazashvili, and P. Coleman, Nature 407, 351 (2000).
  * Prochaska  _et al._ (2020) L. Prochaska, X. Li, D. C. MacFarland, A. M. Andrews, M. Bonta, E. F. Bianco, S. Yazdi, W. Schrenk, H. Detz, A. Limbeck, Q. Si, E. Ringe, G. Strasser, J. Kono, and S. Paschen, Science 367, 285 (2020).
  * Chen  _et al._ (2022) L. Chen, D. T. Lowder, E. Bakali, A. Andrews, W. Schrenk, M. Waas, R. Svagera, G. Eguchi, L. Prochaska, Q. Si, _et al._ , arXiv preprint arXiv:2206.00673 (2022).
  * Keimer and Moore (2017) B. Keimer and J. E. Moore, Nat. Phys. 13, 1045 (2017).
  * Phillips  _et al._ (2022) P. W. Phillips, N. E. Hussey, and P. Abbamonte, arXiv preprint arXiv:2205.12979 (2022).
  * Si  _et al._ (2014) Q. Si, J. H. Pixley, E. Nica, S. J. Yamamoto, P. Goswami, R. Yu, and S. Kirchner, Journal of the Physical Society of Japan 83, 061005 (2014).
  * Si and Smith (1996a) Q. Si and J. L. Smith, Phys. Rev. Lett. 77, 3391 (1996a).
  * Smith and Si (2000) J. L. Smith and Q. Si, Phys. Rev. B 61, 5184 (2000).
  * Chitra and Kotliar (2000) R. Chitra and G. Kotliar, Phys. Rev. Lett. 84, 3678 (2000).
  * Si and Smith (1996b) Q. Si and J. L. Smith, Phys. Rev. Lett. 77, 3391 (1996b).
  * Smith and Si (1999) J. L. Smith and Q. Si, Europhysics Letters (EPL) 45, 228 (1999).
  * Sengupta (2000) A. M. Sengupta, Phys. Rev. B 61, 4041 (2000).
  * Si  _et al._ (2003) Q. Si, S. Rabello, K. Ingersent, and J. L. Smith, Phys. Rev. B 68, 115103 (2003).
  * Zhu and Si (2002) L. Zhu and Q. Si, Phys. Rev. B 66, 024426 (2002).
  * Zaránd and Demler (2002) G. Zaránd and E. Demler, Phys. Rev. B 66, 024427 (2002).
  * Grempel and Si (2003) D. Grempel and Q. Si, Phys. Rev. Lett. 91, 026401 (2003).
  * Zhu  _et al._ (2003) J. Zhu, D. Grempel, and Q. Si, Phys. Rev. Lett. 91, 156404 (2003).
  * Glossop and Ingersent (2007) M. Glossop and K. Ingersent, Phys. Rev. Lett. 99, 227203 (2007).
  * Zhu  _et al._ (2007) J.-X. Zhu, S. Kirchner, R. Bulla, and Q. Si, Phys. Rev. Lett. 99, 227204 (2007).
  * Cai and Si (2019) A. Cai and Q. Si, Phys. Rev. B 100, 014439 (2019).
  * Cuomo  _et al._ (2022) G. Cuomo, Z. Komargodski, M. Mezei, and A. Raviv-Moshe, arXiv preprint arXiv:2202.00040 (2022).
  * (40) Supplementary Materials .
  * Auerbach (2012) A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer Science & Business Media, 2012).
  * Zinn-Justin (2007) J. Zinn-Justin,  _Phase transitions and renormalization group_ (Oxford University Press on Demand, 2007).
  * Nahum (2022) A. Nahum, arXiv preprint arXiv:2202.08431 (2022).
  * Affleck (1995) I. Affleck, arXiv preprint cond-mat/9512099 (1995).
  * Nahum  _et al._ (2015) A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Phys. Rev. X 5, 041048 (2015).
  * Ma and Wang (2020) R. Ma and C. Wang, Phys. Rev. B 102, 020407 (2020).
  * Abanov and Wiegmann (2000) A. Abanov and P. Wiegmann, Nuclear Physics B 570, 685 (2000).
  * Liu  _et al._ (2021) C.-C. Liu, S. Paschen, and Q. Si, arXiv preprint arXiv:2101.01807 (2021).
  * Han  _et al._ (2022) S. Han, D. J. Schultz, and Y. B. Kim, arXiv preprint arXiv:2206.02808 (2022).
  * Fang  _et al._ (2022) Y. Fang, G. Grissonnanche, A. Legros, S. Verret, F. Laliberté, C. Collignon, A. Ataei, M. Dion, J. Zhou, D. Graf, M. J. Lawler, P. A. Goddard, L. Taillefer, and B. J. Ramshaw, Nature Physics 18, 558 (2022).
  * Badoux  _et al._ (2016) S. Badoux, W. Tabis, F. Laliberté, G. Grissonnanche, B. Vignolle, D. Vignolles, J. Béard, D. A. Bonn, W. N. Hardy, R. Liang, N. Doiron-Leyraud, L. Taillefer, and C. Proust, Nature 531, 210 (2016).
  * Chowdhury  _et al._ (2021) D. Chowdhury, A. Georges, O. Parcollet, and S. Sachdev, arXiv preprint arXiv:2109.05037 (2021).

## Supplementary Material

###  .1 A. Renormalization group analysis

We begin by reiterating that our RG procedure is carried out for the
Hamiltonian defined in Eq. 1. The effective actions $S$ and $S_{\chi}$ are
introduced to describe the next-order fluctuations and calculate the
correlations functions.

The relevant correlation functions are:

$\displaystyle\langle\sum_{\mu}[\phi^{\mu}(\bm{r},\tau=0)]^{2}\rangle_{imp}\,,$
$\displaystyle\sum_{\sigma,a}\langle
c_{\sigma,a}(\bm{r_{0}},\tau)c_{\sigma,a}^{\dagger}(\bm{r_{0}},\tau)\rangle_{imp}\,.$
(S1)

Here, $\langle\cdot\rangle_{imp}$ means we subtract the free-bulk
contributions, and the bulk system has spacetime dimension $d=4-\epsilon$.
Before performing calculations, we note that the exact formula of above
correlations functions are

$\displaystyle\langle\sum_{\mu}\phi^{\mu}(\bm{r},\tau=0)^{2}\rangle_{imp}=\int_{\tau,\tau^{\prime}}\frac{A_{\phi}}{2\pi}\frac{h}{(\bm{r}^{2}+\tau^{2})^{(2-\epsilon)/2}(\bm{r}^{2}+\tau^{\prime
2})^{(2-\epsilon)/2}}\langle\sum_{\mu}S^{\mu}(\tau)S^{\mu}(\tau^{\prime})\rangle_{imp}$
$\displaystyle-\sum_{\sigma,a}\langle
c_{\sigma,a}(\bm{r_{0}},\tau)c_{\sigma,a}^{\dagger}(\bm{r_{0}},\tau)\rangle_{imp}=K\bigg{\langle}\bigg{[}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau^{\prime}}+\epsilon_{\bm{k}})+\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}S^{\mu}(\tau)\sigma^{\mu}\bigg{)}\bigg{]}\bigg{\rangle}_{imp}-G_{0}(\bm{r}_{0},\tau),$
(S2)

where

$\displaystyle
h=-4A_{\phi}g^{2}\Gamma(-1+\epsilon)\sin(\frac{\epsilon\pi}{2})\approx[2\pi
A_{\phi}+O(\epsilon)]g^{2}.$ (S3)

$A_{\phi}$ is the normalization factor of boson propagator and
$G_{0}(\bm{r}_{0},\tau)$ is the free fermion propagator. Above equations
suggest the fermion and boson correlation functions can be represented as
functional of $n$-point correlation functions of $S^{\mu}$ fields. Then, it’s
sufficient to evaluate the correlation functions of $S^{\mu}$. To do so, we
first derive the effective action of $S^{\mu}$ by integrating out fermions and
bosons:

$\displaystyle
S_{eff}=-iSS_{B}-K\text{Tr}\log\bigg{(}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau}+\epsilon_{\bm{k}})-\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}S^{\mu}(\tau)\frac{\sigma^{\mu}}{2}\bigg{)}.$
(S4)

Using $S_{eff}$ and Eq. S2, we’re able to calculate fermion and boson
correlation functions. The detailed calculations are shown in Supplementary
Material (SM) B and the results are

$\displaystyle\langle\sum_{\mu}[\phi^{\mu}(\bm{r},\tau=0)]^{2}\rangle_{imp}$
$\displaystyle\propto$
$\displaystyle\frac{1}{|\bm{r}|^{2}}h\bigg{[}1-\frac{1}{S\pi
r}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{S\pi\epsilon}[\frac{\log(1+h/c(J)}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}]\bigg{]}$
$\displaystyle\sum_{\sigma,a}\langle
c_{\sigma,a}(\bm{r_{0}},\tau)c_{\sigma,a}^{\dagger}(\bm{r_{0}},\tau)\rangle_{imp}$
$\displaystyle\propto$
$\displaystyle\frac{J^{2}\text{sgn}(\tau)}{(1+J^{2})|\tau|}\bigg{[}1+\frac{-1-J^{2}}{(1+J^{2})^{2}\pi
Sr}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{(1+J^{2})\pi
Sr}\int_{0}^{J}\frac{2jF(j)}{1+j^{2}}dj$ $\displaystyle-\frac{2}{(1+J^{2})r\pi
S}\int_{0}^{J}\frac{H(j)}{1+j^{2}}dj-\frac{-1+J^{2}}{(1+J^{2})^{2}\pi
S\epsilon}\bigg{(}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}$
$\displaystyle+\frac{2iJ}{(1+J^{2})^{2}\pi\epsilon
S}\bigg{(}-\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}\bigg{]}.$

We are now in the position to carry out renormalization group study. We
introduce the following relation between renormalized couplings $J,h$ and bare
couplings $J_{B},h_{B}$:

$\displaystyle h_{B}=Z_{h}\mu^{\epsilon}h$ $\displaystyle J_{B}=Z_{J}\mu^{r}J$
(S6)

where $Z_{J},Z_{h}$ are coupling renormalization factor, and $\mu$ is a
renormalization energy scale. Without loss of generality, we let

$\displaystyle Z_{h}=1+\frac{1}{S}x_{h}$ $\displaystyle
Z_{J}=1+\frac{1}{S}x_{J}.$

We determine $x_{J}$ and $x_{h}$ by requiring all the poles as a function of
$r,\epsilon$ in Eq. LABEL:eq:corr are minimally removed. This leads to

$\displaystyle x_{h}$ $\displaystyle=$ $\displaystyle\frac{1}{\pi
r}\int_{0}^{J}\frac{F(j)}{j}dj-\frac{1}{\pi\epsilon}\bigg{[}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{]}+O(\frac{1}{S})$
$\displaystyle x_{J}$ $\displaystyle=$ $\displaystyle\frac{1}{2\pi
r}\int_{0}^{J}\frac{F(j)}{j}dj-\frac{1}{\pi
r}\int_{0}^{J}\frac{jF(j)}{1+j^{2}}dj+\frac{1}{r\pi}\int_{0}^{J}\frac{H(j)}{1+j^{2}}dj+\frac{-1+J^{2}}{2(1+J^{2})\pi\epsilon}\bigg{(}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}$
$\displaystyle-\frac{iJ}{(1+J^{2})\pi\epsilon}\bigg{(}-\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}+O(\frac{1}{S})$

Finally, we take $\mu$ derivative of Eq. S6 at fixed $h_{B}$ and $J_{B}$ and
obtain the following beta functions:

$\displaystyle\beta_{J}$ $\displaystyle=$
$\displaystyle\frac{J-J^{3}}{2S\pi(1+J^{2})}F(J)+\frac{J^{2}}{S\pi(1+J^{2})}H(J)+\bigg{[}\frac{-(1-J^{2})Jh+2iJ^{2}h}{2S(1+J^{2})\pi
c(J)(c(J)+h)}+\text{h.c.}\bigg{]}+o(\frac{1}{S})\,,$ $\displaystyle\beta_{h}$
$\displaystyle=$
$\displaystyle-\frac{1}{S}\tilde{\epsilon}h+\frac{h}{S\pi}F(J)-\frac{1}{S\pi}\frac{(c(J)+c(J)^{*})h}{|c(J)|^{2}}-\frac{1}{S\pi}\frac{-2h^{2}-(c(J)+c(J)^{*})h}{(c(J)+h)(c(J)^{*}+h)}+o(\frac{1}{S})\,.$
(S8)

where we also set $r=0$ at the last step.

Note that an alternative way to perform renormalization group study is to
start from the action in Eq. LABEL:eq:action and then expand in powers of
$\chi$ fields Nahum (2022). Combining the quadratic order term $S_{\chi}$ and
fourth-order term $S_{4}$, we could find an interacting field theory of
$\chi$: $S_{\chi}+S_{4}$. Then we are able to perform RG analysis based on
this model. The resulting beta functions are expected to be consistent with
current approach. We reserve such an analysis for a future study.

###  .2 B. Calculations of correlation functions

Here we present the detailed calculation of Eq. S2 using the effective action
$S_{eff}$. We first introduce the spinor representation of $S^{\mu}$ with
$z^{\dagger}\frac{\sigma^{\mu}}{2}z$ and let $z=z_{0}+\delta z$ (see Eq. 3).
This gives us

$\displaystyle S^{\mu}$ $\displaystyle=$ $\displaystyle S\delta_{\mu,z}+\delta
S^{\mu}$ $\displaystyle\delta S^{x}$ $\displaystyle=$
$\displaystyle\frac{i\sqrt{S}}{\sqrt{2}}(\chi^{\dagger}-\chi)+O(\frac{1}{S^{1/2}})$
$\displaystyle\delta S^{y}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{S}}{\sqrt{2}}(\chi^{\dagger}+\chi)+O(\frac{1}{S^{1/2}})$
$\displaystyle\delta S^{z}$ $\displaystyle=$
$\displaystyle-|\chi|^{2}+O(\frac{1}{S}).$ (S9)

We expand the action in Eq. LABEL:eq:action to the quadratic order of $\chi$
fields:

$\displaystyle S_{eff}$ $\displaystyle=$ $\displaystyle S_{0}+S_{\chi}$
$\displaystyle S_{0}$ $\displaystyle=$
$\displaystyle-\frac{Sg^{2}}{4}\int_{\tau,\tau^{\prime}}\frac{\sum_{\mu}(z_{0}^{\dagger}\sigma^{\mu}z_{0})^{2}}{2}G_{\phi}(\tau-\tau^{\prime})-K\text{Tr}\log\bigg{(}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau^{\prime}}+\epsilon_{\bm{k}})+\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}(z_{0}^{\dagger}\sigma^{\mu}z_{0})\sigma^{\mu}\bigg{)},$
(S10)

with $S_{0}$ the saddle point contribution at order of $S$ and $S_{\chi}$
describing fluctuations at order $S^{0}$. We now derive the formula of
$S_{\chi}$. It has contribution from three part: the Berry phase, bosonic
coupling and Kondo interaction. The Berry phase term gives

$\displaystyle\int i\Omega|\chi(\Omega)|^{2}\frac{d\Omega}{2\pi}\,,$ (S11)

while the bosonic coupling yields

$\displaystyle-\frac{g^{2}}{2}\int_{\tau,\tau^{\prime}}G_{\phi}(\tau-\tau^{\prime})[\chi^{\dagger}(\tau)\chi(\tau^{\prime})+\chi(\tau)\chi^{\dagger}(\tau^{\prime})-|\chi(\tau)|^{2}-|\chi(\tau^{\prime})|^{2}]$
(S12) $\displaystyle=$ $\displaystyle\int
h|\Omega|^{1-\epsilon}|\chi(i\Omega)|^{2}\frac{d\Omega}{2\pi}\,,$

and, finally, the Kondo coupling gives (by expanding the fermion determinant)

$\displaystyle\frac{\kappa
J_{K}^{2}}{4}\int_{\tau,\tau^{\prime}}\bigg{[}G_{s}(\tau-\tau^{\prime})G_{s}(\tau^{\prime}-\tau)-G_{a}(\tau-\tau^{\prime})G_{a}(\tau^{\prime}-\tau)\bigg{]}\bigg{[}\chi^{\dagger}(\tau)\chi(\tau^{\prime})+\chi^{\dagger}(\tau^{\prime})\chi(\tau)-|\chi(\tau)|^{2}-|\chi(\tau^{\prime})|^{2}\bigg{]}$
(S13) $\displaystyle+\frac{\kappa
J_{K}^{2}}{2}\int_{\tau,\tau^{\prime}}\bigg{[}G_{s}(\tau-\tau^{\prime})G_{a}(\tau^{\prime}-\tau)\bigg{]}\bigg{[}-\chi^{\dagger}(\tau)\chi(\tau^{\prime})+\chi(\tau)\chi^{\dagger}(\tau^{\prime})\bigg{]}$
$\displaystyle=$ $\displaystyle\int\frac{2\kappa
J^{2}|\Omega|^{1-2r}}{\pi(1+J^{2}|\Omega|^{-2r})^{2}}|\chi(i\Omega)|^{2}\frac{d\Omega}{2\pi}+\int\frac{i\Omega\kappa}{\pi}\bigg{[}\frac{J|\Omega|^{-r}-J^{3}|\Omega|^{-3r}}{(1+J^{2}|\Omega|^{-2r})^{2}}-\arctan(J|\Omega|^{-r})\bigg{]}|\chi(i\Omega)|^{2}\frac{d\Omega}{2\pi},$

where $G_{s}(\tau)=\sum_{\sigma}G_{\sigma}(\tau)$ and
$G_{a}(\tau)=\sum_{\sigma}\sigma G_{\sigma}(\tau)$ and $G_{\sigma}(\tau)$
denotes the conduction electron Green’s function at leading order $(S^{1})$.
The explicit formula of these functions are

$\displaystyle G_{s}(i\omega)$
$\displaystyle=\frac{-2i\rho_{0}\pi|\omega|^{-r}}{1+\rho_{0}^{2}\pi^{2}J_{K}^{2}|\omega|^{-2r}}\text{sgn}(\omega)$
$\displaystyle G_{a}(i\omega)$
$\displaystyle=\frac{-2\rho_{0}^{2}\pi^{2}J_{K}|\omega|^{-2r}}{1+\rho_{0}^{2}\pi^{2}J_{K}^{2}|\omega|^{-2r}}.$
(S14)

In addition, the coefficient before each $|\Omega|^{-r}$ should be a function
of $r$. However, we only keep the leading order term ($r^{0}$) since we would
set $r=0$ in the final step. Combining Eqs. S11, S12 and S13, we have the
quadratic theory of the $\chi$ fields:

$\displaystyle S_{\chi}=\int
G^{-1}_{\chi}(i\Omega)|\chi(i\Omega)|^{2}\frac{d\Omega}{2\pi}$ $\displaystyle
G_{\chi}^{-1}(i\Omega)=i\Omega\bigg{[}1+p(J|\Omega|^{-r})\bigg{]}+|\Omega|\bigg{[}q(J|\Omega|^{-r})+h|\Omega|^{-\epsilon}\bigg{]}.$
(S15)

For future reference, we introduce the series expansion of $G_{\chi}(i\Omega)$

$\displaystyle G_{\chi}(i\Omega)$ $\displaystyle=$
$\displaystyle\frac{1}{|\Omega|}\frac{1}{i(\text{sgn}(\Omega)+p(J|\Omega|^{-r}))+q(J|\Omega|^{-r})}+\frac{1}{|\Omega|}\sum_{m=1}^{\infty}\frac{(-h)^{m}|\Omega|^{-m\epsilon}}{[i(\text{sgn}(\Omega)+p(J))+q(J)]^{m+1}}$
(S16)

Since we set $r=0$ in the last step, we drop the $r$ dependence in the
summation by setting $r=0$. Including this contribution won’t change the beta
functions at the leading order. We can further expand the real and imaginary
part of $G_{\chi}(i\Omega)$ as the following

$\displaystyle
G_{\chi}(i\Omega)+G_{\chi}(-i\Omega)=\sum_{n=0}^{\infty}a_{n}J^{n}|\Omega|^{-1-nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}$
$\displaystyle
G_{\chi}(i\Omega)-G_{\chi}(i\Omega)=\sum_{n=0}^{\infty}ib_{n}J^{n}|\Omega|^{1-nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{1-m\epsilon}$
(S17)

where $a_{n},b_{n}$ are the Taylor coefficients of $F(j)$ and $H(j)$ at $j=0$
and

$\displaystyle F(j)$ $\displaystyle=$
$\displaystyle\frac{2q(j)}{(1+p(j))^{2}+q(j)^{2}}$ $\displaystyle H(j)$
$\displaystyle=$ $\displaystyle\frac{-2(1+p(j))}{(1+p(j))^{2}+q(j)^{2}}$
$\displaystyle c(J)$ $\displaystyle=$ $\displaystyle i(1+p(J))+q(J).$

Now, we combine Eqs. S2, S9, S10 and S15 and calculate the correlation
functions of fermions and bosons. We first expand Eq. S2 using Eq. S9. For
bosons, we have

$\displaystyle\langle\sum_{\mu}\phi^{\mu}(\bm{r},\tau=0)^{2}\rangle_{imp}$
$\displaystyle=$
$\displaystyle\int_{\tau,\tau^{\prime}}\frac{A_{\phi}}{2\pi}\frac{h}{(\bm{r}^{2}+\tau^{2})^{(2-\epsilon)/2}(\bm{r}^{2}+\tau^{\prime
2})^{(2-\epsilon)/2}}S^{2}$
$\displaystyle+\int_{\tau,\tau^{\prime}}\frac{A_{\phi}}{2\pi}\frac{h}{(\bm{r}^{2}+\tau^{2})^{(2-\epsilon)/2}(\bm{r}^{2}+\tau^{\prime
2})^{(2-\epsilon)/2}}\langle-S\delta S^{z}(\tau)-S\delta
S^{z}(\tau^{\prime})+\delta S^{x}(\tau)\delta S^{x}(\tau^{\prime})+\delta
S^{y}(\tau)\delta S^{y}(\tau^{\prime})\rangle_{S_{\chi}}+O(S^{0})$

where the last line denotes the contribution of the next-order fluctuations
and are evaluated with respect to $S_{\chi}$. For the fermionic fields, it
gives

$\displaystyle-\langle\sum_{a,\sigma}c_{\bm{r_{0}},a,\sigma}(\tau)c_{\bm{r_{0}},a,\sigma}^{\dagger}(0)\rangle$
(S19) $\displaystyle=$ $\displaystyle K\bigg{[}G_{s}(\tau)-G_{0}(\tau)$
$\displaystyle+\int_{\tau_{1}}\sum_{\mu}\frac{J_{K}\text{Tr}[G(\tau-\tau_{1})\sigma^{\mu}G(\tau_{1})]}{S}\langle\delta
S^{\mu}(\tau_{1})\rangle_{S_{\chi}}+\int_{\tau_{1},\tau_{2}}\sum_{\mu,\nu}\frac{J^{2}_{K}\text{Tr}[G(\tau-\tau_{1})\sigma^{\mu}G(\tau_{1}-\tau_{2})\sigma^{\nu}G(\tau_{2})]}{S^{2}}\langle\delta
S^{\mu}(\tau_{1})\delta S^{\nu}(\tau_{2})\rangle_{S_{\chi}}\bigg{]}$
$\displaystyle+O(S^{0})$

where the third line denotes the contribution of the next-order fluctuations
and are evaluated with respect to $S_{\chi}$ and
$G(\tau)=\text{diag}[\frac{G_{s}(\tau)+G_{a}(\tau)}{2},\frac{G_{s}(\tau)-G_{a}(\tau)}{2}]$.

To get the final results, we only need to calculate the correlation functions
of $\delta S^{\mu}$ (or $\chi$) that appearing in the above equations. The
boson correlation function becomes

$\displaystyle\langle\sum_{\mu}\phi^{\mu}(\bm{r},\tau=0)^{2}\rangle_{imp}$
(S20) $\displaystyle=$
$\displaystyle\int_{\tau,\tau^{\prime}}\frac{A_{\phi}}{2\pi}\frac{h}{(\bm{r}^{2}+\tau^{2})^{(2-\epsilon)/2}(\bm{r}^{2}+\tau^{\prime
2})^{(2-\epsilon)/2}}S^{2}$
$\displaystyle\bigg{[}1+\frac{1}{S}\langle\chi^{\dagger}(\tau)\chi(\tau^{\prime})+\chi^{\dagger}(\tau^{\prime})\chi(\tau^{\prime})-|\chi(\tau)|^{2}-|\chi(\tau^{\prime})|^{2}\rangle\bigg{]}+O(S^{0})$
$\displaystyle=$ $\displaystyle\frac{A_{\phi}\pi
S^{2}}{2|\bm{r}|^{2}}h+\frac{A_{\phi}Sh}{2\pi}\int_{0}^{\infty}2(G_{\chi}(\Omega)+G_{\chi}(-\Omega))\bigg{(}2\pi
x^{-(1-\epsilon)}B[\frac{\epsilon-1}{2},x]^{2}|\Omega|^{2-2\epsilon}-\frac{\pi^{2}}{|\bm{r}|^{2}}\bigg{)}\frac{d\Omega}{2\pi}+O(S^{0})$
$\displaystyle=$ $\displaystyle\frac{A_{\phi}\pi
S^{2}}{2|\bm{r}|^{2}}h+\frac{A_{\phi}Sh}{2\pi}\int_{0}^{\infty}2\bigg{[}\sum_{n=0}^{\infty}a_{n}J^{n}|\Omega|^{1-nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}\bigg{]}$
$\displaystyle\frac{1}{|\bm{r}|^{2}}\bigg{(}2\pi
x^{1-\epsilon}B[\frac{\epsilon-1}{2},x]^{2}-\pi
2^{-\epsilon}\Gamma((1-\epsilon)/2)^{2}\bigg{)}\frac{d\Omega}{2\pi}+O(S^{0})$
$\displaystyle=$ $\displaystyle\frac{A_{\phi}\pi
S^{2}}{2|\bm{r}|^{2}}h+\frac{A_{\phi}Sh}{\pi|\bm{r}|^{2}}\bigg{[}\frac{1}{r}\sum_{n=0}^{\infty}\frac{-a_{n}J^{n}\pi}{2n}+\frac{1}{\epsilon}\sum_{m=1}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}\frac{-\pi}{2n\epsilon}+O(r^{0},\epsilon^{0})\bigg{]}+O(S^{0})$
$\displaystyle=$ $\displaystyle\frac{A_{\phi}\pi
S^{2}}{2|\bm{r}|^{2}}h\bigg{\\{}1-\frac{1}{S}\bigg{[}\frac{1}{\pi
r}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{\pi\epsilon}[\frac{\log(1+h/c(J)}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}]+O(\epsilon^{0},r^{0})\bigg{]}+o(\frac{1}{S})\bigg{\\}}$

where $x=|\Omega||\bm{r}|$ and $B(a,x)$ is the modified Bessel functions of
the second kind. In the above expression, we only keep track of the poles as a
function of $\epsilon$ and $r$. We also drop the $\frac{r}{\epsilon}$ poles
which goes to zero as we set $r=0$.

For the fermion correlation function, we first perform Fourier transformation
and separate the next-order contributions into four parts
$I_{1},I_{2},I_{3},I_{4}$:

$\displaystyle
G_{f}(i\omega)=-\int_{\tau}\langle\sum_{a,\sigma}c_{\bm{r_{0}},a,\sigma}(\tau)c_{\bm{r_{0}},a,\sigma}^{\dagger}(0)\rangle_{imp}e^{i\omega\tau}=K\bigg{[}\frac{2i\pi\rho_{0}J^{2}}{1+J^{2}}\text{sgn}(\omega)+I_{1}+I_{2}+I_{3}+I_{4}+o(\frac{1}{S})\bigg{]}$
(S21)

Each part is shown below:

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{J_{K}^{2}}{4S}(G_{s}(i\omega)^{2}+G_{a}(i\omega)^{2})\int
G_{\chi}(i\Omega)(G_{s}(i\omega+i\Omega)+G_{s}(i\omega-i\Omega)-G_{s}(i\omega)-G_{s}(i\omega))\frac{d\Omega}{2\pi}$
$\displaystyle=$ $\displaystyle\frac{-2i\rho_{0}\pi
J^{2}(-1+J^{2})}{(1+J^{2})^{2}S}\int_{0}^{\infty}\bigg{[}\sum_{n=0}^{\infty}a_{n}J^{n}|\Omega|^{-1-nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}\bigg{]}$
$\displaystyle\sum_{p=0}^{\infty}(-J^{2})^{p}\bigg{[}\text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r}+\text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}-2\text{sgn}(\omega)|\omega|^{-2pr-r}\bigg{]}\frac{d\Omega}{2\pi}+O(r^{0})$
$\displaystyle=$ $\displaystyle\frac{-i\rho_{0}\pi
J^{2}(-1+J^{2})}{(1+J^{2})^{2}\pi
S}\bigg{[}\sum_{n=0}^{\infty}a_{n}J^{n}\sum_{p=0}^{\infty}(-J^{2})^{P}\frac{-2}{nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}\frac{-2}{m\epsilon}\bigg{]}\text{sgn}(\omega)$
$\displaystyle=$ $\displaystyle\frac{2i\rho_{0}\pi
J^{2}(-1+J^{2})}{(1+J^{2})^{3}\pi
S}\bigg{[}\frac{1}{r}\int_{0}^{J}\frac{F(j)}{j}dj-\frac{1}{\epsilon}\bigg{(}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}+O(r^{0},\epsilon^{0})\bigg{)}\bigg{]}\text{sgn}(\omega).$

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle-\frac{J_{K}^{2}}{2S}G_{s}(i\omega)G_{a}(i\omega)\int
G_{\chi}(i\Omega)(G_{a}(i\omega+i\Omega)+G_{a}(i\omega-i\Omega)-G_{a}(i\omega)-G_{a}(i\omega))\frac{d\Omega}{2\pi}$
$\displaystyle=$
$\displaystyle-\frac{2iJ^{3}(-2J\rho_{0}\pi)}{(1+J^{2})^{2}S}\text{sgn}(\omega)\int_{0}^{\infty}\bigg{[}\sum_{n=0}^{\infty}a_{n}J^{n}|\Omega|^{-1-nr}+\sum_{m=1}^{\infty}(-h)^{m}\bigg{(}\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}\bigg{]}$
$\displaystyle\sum_{p=0}^{\infty}(-J^{2})^{p}\bigg{[}|\omega+\Omega|^{-2pr-2r}+|\omega-\Omega|^{-2pr-2r}-2|\omega|^{-2pr-2r}\bigg{]}\frac{d\Omega}{2\pi}$
$\displaystyle=$
$\displaystyle-\frac{2iJ^{3}(-2J\rho_{0}\pi)}{2\pi(1+J^{2})^{2}S}\text{sgn}(\omega)\bigg{[}\sum_{n=0}^{\infty}a_{n}\sum_{p=0}^{\infty}(-J^{2})^{p}J^{n}(\frac{2}{2pr+2r+nr}-\frac{2}{nr})+O(r^{0},\epsilon^{0})\bigg{]}$
$\displaystyle=$
$\displaystyle\frac{2iJ^{4}\rho_{0}\pi}{\pi(1+J^{2})^{2}S}\text{sgn}(\omega)\bigg{[}\frac{1}{rJ^{2}}\int_{0}^{J}\frac{2jF(j)}{1+j^{2}}dj-\frac{2}{r(1+J^{2})}\int_{0}^{J}\frac{F(j)}{j}dj+O(r^{0},\epsilon^{0})\bigg{]}.$

$\displaystyle I_{3}$ $\displaystyle=$
$\displaystyle\frac{J_{K}^{2}}{4S}\bigg{[}G_{a}(i\omega)^{2}+G_{s}(i\omega)^{2}\bigg{]}\int
G_{\chi}(i\Omega)[G_{a}(i\Omega+i\omega)-G_{a}(i\omega-i\Omega)]\frac{d\Omega}{2\pi}$
$\displaystyle=$
$\displaystyle\frac{J^{2}(-1+J^{2})}{2\pi(1+J^{2})^{2}S}\int_{0}^{\infty}\bigg{[}\sum_{n=0}ib_{n}J^{n}|\Omega|^{1-nr}+\sum_{m=1}(-h)^{m}\bigg{(}-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}\bigg{]}$
$\displaystyle\bigg{[}-2J\rho_{0}\pi\sum_{p=0}^{\infty}(-J^{2})^{p}\bigg{(}|\Omega+\omega|^{-2r-2pr}-|-\Omega+\omega|^{-2r-2pr}\bigg{)}\bigg{]}d\Omega$
$\displaystyle=$ $\displaystyle\frac{1}{S}[0+O(r^{0},\epsilon^{0})].$

$\displaystyle I_{4}$ $\displaystyle=$
$\displaystyle\frac{-J_{K}^{2}}{2}G_{a}(i\omega)G_{s}(i\omega)\int
G_{\chi}(i\Omega)\bigg{(}G_{s}(i\Omega+i\omega)-G_{s}(i\omega-i\Omega)\bigg{)}\frac{d\Omega}{2\pi}$
$\displaystyle=$
$\displaystyle\frac{-2J^{3}\rho_{0}\pi}{S(1+J^{2})^{2}}\int_{0}^{\infty}\bigg{[}\sum_{n=0}ib_{n}J^{n}|\Omega|^{1-nr}+\sum_{m=1}(-h)^{m}\bigg{(}-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}|\Omega|^{-1-m\epsilon}\bigg{]}$
$\displaystyle\sum_{p=0}(-J^{2})^{p}\bigg{[}\text{sgn}(\omega+\Omega)|\omega+\Omega|^{-2pr-r}-\text{sgn}(\omega-\Omega)|\omega-\Omega|^{-2pr-r}\bigg{]}\frac{d\Omega}{\pi}$
$\displaystyle=$
$\displaystyle\frac{-2J^{3}\rho_{0}\pi}{S(1+J^{2})^{2}}\bigg{[}\sum_{n=0}ib_{n}J^{n}\sum_{p=0}(-J^{2})^{p}\frac{2}{\pi(nr+2pr+r)}+\sum_{m=1}(-h)^{m}\frac{1}{1+J^{2}}\bigg{(}-\frac{1}{(c(J)^{*})^{m+1}}+\frac{1}{c(J)^{m+1}}\bigg{)}\frac{2}{\pi
m\epsilon}\bigg{]}\text{sgn}(\omega)$ $\displaystyle=$
$\displaystyle\frac{-2J^{3}\rho_{0}\pi}{(1+J^{2})^{2}}\bigg{[}\frac{2}{r\pi
J}\int_{0}^{J}\frac{iH(j)}{(1+j^{2})}dj+\frac{2}{\pi\epsilon(1+J^{2})}\bigg{(}-\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}+O(r^{0},\epsilon^{0})\bigg{]}\text{sgn}(\omega).$

In the above derivation, we only keep terms that diverge as $r\rightarrow 0$
and $\epsilon\rightarrow 0$, with $r$ going to zero first.

Summing over all the contributions, we have

$\displaystyle G_{f}(i\omega)$ $\displaystyle=$ $\displaystyle
2iK\pi\rho_{0}\text{sgn}(\omega)\frac{J^{2}}{1+J^{2}}\bigg{[}1+\frac{-1+J^{2}}{(1+J^{2})^{2}\pi
Sr}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{(1+J^{2})\pi
Sr}\int_{0}^{J}\frac{2jF(j)}{1+j^{2}}dj$
$\displaystyle-\frac{2J^{2}}{(1+J^{2})^{2}\pi
Sr}\int_{0}^{J}\frac{F(j)}{j}dj-\frac{2}{(1+J^{2})r\pi
S}\int_{0}^{J}\frac{H(j)}{1+j^{2}}dj$
$\displaystyle-\frac{-1+J^{2}}{(1+J^{2})^{2}\pi
S\epsilon}\bigg{(}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}+\frac{2iJ}{(1+J^{2})^{2}\pi\epsilon
S}\bigg{(}-\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}\bigg{]}$

In summary, we have

$\displaystyle\langle\sum_{\mu}[\phi^{\mu}(\bm{r},\tau=0)]^{2}\rangle_{imp}$
$\displaystyle\propto$
$\displaystyle\frac{1}{|\bm{r}|^{2}}h\bigg{[}1-\frac{1}{S\pi
r}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{S\pi\epsilon}[\frac{\log(1+h/c(J)}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}]\bigg{]}$
$\displaystyle\sum_{\sigma,a}\langle
c_{\sigma,a}(\bm{r_{0}},\tau)c_{\sigma,a}^{\dagger}(\bm{r_{0}},\tau)\rangle_{imp}$
$\displaystyle\propto$
$\displaystyle\frac{J^{2}\text{sgn}(\tau)}{(1+J^{2})|\tau|}\bigg{[}1+\frac{-1-J^{2}}{(1+J^{2})^{2}\pi
Sr}\int_{0}^{J}\frac{F(j)}{j}dj+\frac{1}{(1+J^{2})\pi
Sr}\int_{0}^{J}\frac{2jF(j)}{1+j^{2}}dj$ $\displaystyle-\frac{2}{(1+J^{2})r\pi
S}\int_{0}^{J}\frac{H(j)}{1+j^{2}}dj-\frac{-1+J^{2}}{(1+J^{2})^{2}\pi
S\epsilon}\bigg{(}\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}$
$\displaystyle+\frac{2iJ}{(1+J^{2})^{2}\pi\epsilon
S}\bigg{(}-\frac{\log(1+h/c(J))}{c(J)}+\frac{\log(1+h/c(J)^{*})}{c(J)^{*}}\bigg{)}\bigg{]}.$

###  .3 C. Beta function analysis

We numerically evaluate the beta functions and the resulting RG flows are
shown in Fig. S1. Besides numerical evaluations, we analyze the beta equations
in certain extreme cases.

At $h=0$, the model becomes a standard Kondo impurity model. The beta function
of $J$ is

$\displaystyle\beta_{J}\bigg{|}_{h=0}$ $\displaystyle=$
$\displaystyle\frac{J-J^{3}}{2S\pi(1+J^{2})}F(J)+\frac{J^{2}}{S\pi(1+J^{2})}H(J)$
(S23) $\displaystyle=$ $\displaystyle-\frac{2J^{2}}{S\pi}+\frac{2\kappa
J^{3}}{S\pi^{2}}+O(J^{4}/S).$

Numerically, we find $J$ is always relevant at $\kappa=2$ and the system flows
to a strong coupling fixed point $K$ at $J=\infty$. The first two terms in the
second line are also consistent with the result obtained from perturbation
theory as shown in Ref. Affleck (1995) (note that we’ve scaled the Kondo
coupling).

At $J=0$, the model becomes a Bose-Kondo model with the following beta
function of $h$:

$\displaystyle\beta_{h}\bigg{|}_{J=0}=-\frac{1}{S}\tilde{\epsilon}h+\frac{2h^{2}}{S\pi(1+h^{2})}.$
(S24)

This is consistent with previous results obtained in Refs. Cuomo  _et al._
(2022); Nahum (2022). It gives two stable fixed points at
$h=(1-\sqrt{1-\tilde{\epsilon}^{2}})/{\tilde{\epsilon}},\infty$ (L and L′) and
one critical point at $h=(1+\sqrt{1-\tilde{\epsilon}^{2}})/{\tilde{\epsilon}}$
(LC). After introducing the Kondo coupling, we find $J$ is irrelevant at both
L and LC, but marginal at L′.

At $h=\infty$, we have a line of fixed points (L′) where $J$ is always
marginal. The beta function at the large $h$ is:

$\displaystyle\beta_{J}=\frac{J-J^{3}}{\pi
S(1+J^{2})}\frac{1}{h}+O(\frac{1}{h^{2}S}).$ (S25)

Clearly, $\beta_{J}\rightarrow 0$ as $h\rightarrow\infty$.

At $J=\infty$, besides the K fixed point, there is an unstable fixed point KD′
at $(h_{c}=\frac{2}{\tilde{\epsilon}\pi},J=\infty)$. To see the existence of
such a critical point, we give the beta function of $h$ at large $J$:

$\displaystyle\beta_{h}=-\frac{\tilde{\epsilon}h}{S}+\frac{2}{S\pi}+O(\frac{1}{SJ}).$
(S26)

It suggests that $h$ flows to $0$ ($\infty$) when $h<(>)h_{c}$ at
$J\rightarrow\infty$.

Figure S1: Renormalization group flow at
$\tilde{\epsilon}<\tilde{\epsilon}_{1}$ (a),
$\tilde{\epsilon}_{1}<\tilde{\epsilon}<\tilde{\epsilon}_{2}$ (b) and
$\tilde{\epsilon}_{2}<\tilde{\epsilon}$ (c) obtained from numerical
evaluation. For purpose of illustration, we normalize all the vectors to have
the same length.

###  .4 D. Topological contribution from the Kondo coupling

In this section, we’ll show that the Kondo coupling contributes an additional
topological term to the impurity spin at sufficiently large values of $J_{K}$.
For the purpose of discussion, we consider pure Kondo model, in other words,
the BFK model at $g=0$ or $h=0$. We then integrate out the fermionic fields,
which gives the following action

$\displaystyle
S_{K}=-iSS_{B}-K\text{Tr}\log\bigg{(}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau}+\epsilon_{\bm{k}})-\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}S^{\mu}(\tau)\frac{\sigma^{\mu}}{2}\bigg{)}.$
(S27)

We now evaluate the fermion determinant contribution at sufficiently large
$J_{K}$ via a Trotter decomposition

$\displaystyle\exp\bigg{(}K\text{Tr}\log\bigg{(}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau}+\epsilon_{\bm{k}})-\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}S^{\mu}(\tau)\frac{\sigma^{\mu}}{2}\bigg{)}\bigg{)}=\bigg{[}\lim_{N\rightarrow\infty}\text{Tr}[\prod^{N}_{i}\exp[-H(\tau_{i})d\tau]]\bigg{]}^{K}$
(S28)

where

$\displaystyle
H(\tau_{i})=J_{K}\sum_{\mu}c_{\bm{r}_{0}}\frac{\sigma^{\mu}}{2}c_{\bm{r}_{0}}S^{\mu}(\tau_{i})+\sum_{\bm{k}\sigma}\epsilon_{\bm{k}}c_{\bm{k},\sigma}^{\dagger}c_{\bm{k},\sigma}.$
(S29)

Since the $K$ orbitals are degenerate, it’s sufficient to consider one orbital
and multiply the contribution by $K$. When $J_{K}$ is sufficiently large, the
hopping between site $\bm{r}_{0}$ and the other sites is a high order effect
($\frac{t^{2}}{J_{K}}$). Here, we focus on the leading order contribution, and
it’s sufficient to take
$H(\tau_{i})=J_{K}/2\sum_{\mu}c_{\bm{r}_{0}}\sigma^{\mu}c_{\bm{r}_{0}}S^{\mu}(\tau_{i})$.
The eigenstates and eigenvalues of $H(\tau_{i})$ are

$\displaystyle|v_{+}(\tau_{i})\rangle=\cos(\frac{\theta(\tau_{i})}{2})c_{\bm{r_{0}},\uparrow}^{\dagger}|0\rangle+\sin(\frac{\theta(\tau_{i})}{2})e^{i\phi(\tau_{i})}c_{\bm{r_{0}},\downarrow}^{\dagger}|0\rangle$
$\displaystyle|v_{-}(\tau_{i})\rangle=-\sin(\frac{\theta(\tau_{i})}{2})c_{\bm{r_{0}},\uparrow}^{\dagger}|0\rangle+\cos(\frac{\theta(\tau_{i})}{2})e^{i\phi(\tau_{i})}c_{\bm{r_{0}},\downarrow}^{\dagger}|0\rangle$
$\displaystyle|v_{0}(\tau_{i})\rangle=|0\rangle$
$\displaystyle|v_{d}(\tau_{i})\rangle=|\uparrow\downarrow\rangle$
$\displaystyle E_{\pm}(\tau_{i})=\pm\frac{J_{K}S}{2},\quad
E_{0}(\tau_{i})=0,\quad E_{d}(\tau_{i})=0$ (S30)

where we let
$\bm{S}(\tau_{i})=S(\sin(\theta(\tau_{i}))\cos(\phi(\tau_{i})),\sin(\theta(\tau_{i}))\sin(\phi(\tau_{i})),\cos(\theta(\tau_{i})))$.
We have the following decomposition via the eigenvalues and eigenstates

$\displaystyle e^{-H(\tau_{i})d\tau}$ $\displaystyle=$ $\displaystyle
e^{-E_{+}(\tau_{i})d\tau}|v_{+}(\tau_{i})\rangle\langle
v_{+}(\tau_{i})|+e^{-E_{-}(\tau_{i})d\tau}|v_{-}(\tau_{i})\rangle\langle
v_{-}(\tau_{i})|$ (S31)
$\displaystyle+e^{-E_{0}(\tau_{i})d\tau}|v_{0}(\tau_{i})\rangle\langle
v_{0}(\tau_{i})|+e^{-E_{d}(\tau_{i})d\tau}|v_{d}(\tau_{i})\rangle\langle
v_{d}(\tau_{i})|.$

Combining the above identity with a Trotter decomposition of the fermion
determinant, we obtain the following leading order contribution

$\displaystyle\exp(K\text{Tr}\log\bigg{(}\delta_{\tau,\tau^{\prime}}\delta_{k,k^{\prime}}(\partial_{\tau}+\epsilon_{\bm{k}})-\delta_{\tau,\tau^{\prime}}J_{K}\sum_{\mu}S^{\mu}(\tau)\frac{\sigma^{\mu}}{2}\bigg{)})$
(S32) $\displaystyle\propto$
$\displaystyle\exp\bigg{(}K\int_{\tau}\langle\partial_{\tau}v_{-}(\tau)|v_{-}(\tau)\rangle
d\tau+O(\frac{1}{J_{K}})\bigg{)}$ $\displaystyle\propto$
$\displaystyle\exp\bigg{(}i\frac{K}{2}\int_{\tau}(1-\cos(\theta))\partial_{\tau}\phi
d\tau+O(\frac{1}{J_{K}})\bigg{)}$ $\displaystyle=$
$\displaystyle\exp\bigg{(}-i\frac{K}{2}S_{B}+O(\frac{1}{J_{K}})\bigg{)}.$

The effective action then becomes

$\displaystyle S_{K}=-i(S-\frac{K}{2})S_{B}+O(\frac{1}{J_{K}})$ (S33)

which indicates the non-trivial topological contribution from the coupling to
the gapless fermions. Qualitatively, the spin has been Kondo-screened, with
the effective spin becoming $S-\frac{K}{2}$. For the perfect screening case
that we consider, $S-\frac{K}{2}=0$, and we obtain a Kondo singlet.


# Line Spectral Estimation via Unlimited Sampling

###### Abstract

Line spectral estimation (LSE) is a fundamental problem in signal processing
due to its wide applications. For signals that are orders of magnitude larger
than the dynamic range of the analog-to-digital (ADC) threshold, conventional
ADC will clip or saturate, leading to significant information loss. The
Unlimited Sensing Framework (USF) was introduced to avoid saturation through
sampling of the signal modulo. Motivated by the USF, we study the LSE from
modulo samples. By exploiting oversampling and the bounded spectral property,
the US-LSE is proposed to recover the folding instants and perform LSE. Our
numerical simulations show that for oversampling factor $\gamma\geq 10$, the
US-LSE is more stable in a lower signal to noise scenario, in the range of
$20\sim 30$ dB, compared to the existing algorithm. Besides, we process the
real data generated by AWR1642, and show that US-LSE estimates the ranges of
two corners with SNRs of $12$ dB and $23$ dB for oversampling factor
$\gamma=10$, and the normalized dynamic range ${\rm
DR}=\beta_{g}/\lambda\approx 3$, where $\beta_{g}$ is the infinity norm of the
signal.

keywords: Unlimited sampling, orthogonal matching pursuit, line spectral
estimation, modulo samples

##  1 Introduction

Line spectral estimation (LSE) aims to estimate a superposition of several
complex exponentials, which arises in several areas of science and engineering
fields. Many methods such as the fast Fourier transform (FFT), the subspace
based method and the compressed sensing based method have been proposed to do
so [1, 2, 3]. These techniques rely on the assumption that the signals’
dynamic range is within that of the analog-to-digital converters (ADCs). In
practice, however, a signal’s amplitude may be unknown and much larger than
the dynamic range of the ADC. This results in clipping and leads to a new form
of information loss.

In [4], the Unlimited Sampling Framework (USF) was proposed for signal sensing
and recovery. By applying a modulo operation to the signal before the ADC, the
saturation problem was avoided. However, the amplitudes were folded with
values in the range $[-\lambda,\lambda]$, where $\lambda>0$ is the ADC
threshold. It has been proven that when the sampling rate is just slightly
above the Nyquist rate, a finite energy bandlimited function is identifiable
by its modulo samples. In addition, a sufficient condition in which the
oversampling rate is greater than $2\pi e$ was provided for recovery. The
proposed recovery algorithm is sensitive to noise due to the use of higher-
order differences. The USF was then further extended to different problems and
signal models such as LSE, DOA estimation, wavelets, mmWave FMCW radar,
finite-rate-of-innovation signals, multi-dimensional signals, quantized
measurements, sparse vector recovery, computed tomography, and graph signals,
etc.[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. In [21], a
periodic function reconstruction was studied. Instead of using higher-order
differences in the time domain for recovery [4], a first-order difference was
adopted and a new Fourier domain recovery algorithm was proposed. This
approach makes the reconstruction more stable and its performance was
validated both in theory and via extensive experiments on their prototype ADC.
In [22], exploiting the fact that bandlimited signals are time-limited, a
robust algorithm was proposed to recover the signal from the modulo samples in
the frequency domain. It was numerically shown that the proposed method has
lower error compared to existing approaches for a given sampling rate, noise
level, and dynamic range of the ADC.

LSE in the noiseless scenario was studied in [5]. Utilizing the general
recovery method [4] for bandlimited signals, the original signal was
recovered, and exponential fitting methods were used to perform LSE. In
contrast, we utilize the fact that both a first-order difference of the simple
function [21] and the line spectra in the frequency domain are sparse. By
exploiting oversampling, the first-order difference of the modulo samples has
spectral beyond bandwidth, which we used to recover the first-order difference
of the simple function via the orthogonal matching pursuit (OMP).

In this work, we address the LSE estimation via USF. By restricting the
frequency on the grid, the measurement model relating to the amplitudes of
frequencies, the modulo samples and the simple function is established. OMP is
then employed to recover the first-order difference of the simple function,
which is further used to recover the spectrum. Numerical simulations involving
real data acquired by mmWave FMCW radar are conducted to demonstrate the
performance of the proposed algorithm, compared with the existing approach.

##  2 Problem Setup and Recovery Algorithm

Consider the continuous line spectral estimation problem

$\displaystyle g(t)=x(t)+w(t),$ (1)

where $x(t)=\sum\limits_{k=1}^{K}a_{k}{\rm e}^{{\rm j}2\pi f_{k}t}$, $f_{k}$
denotes the frequency of the $k$th signal and $0\leq f_{k}\leq f_{\rm max}$,
$t\in{\mathbb{R}}$, $w(t)$ is the noise. Obviously, $f_{\rm max}$ is the
Nyquist frequency, i.e., $f_{\rm Nyq}=f_{\rm max}$. Let $\gamma$ denote the
oversampling factor with respect to the Nyquist frequency $f_{\rm Nyq}$. The
sampling frequency $f_{\rm s}$ is

$\displaystyle f_{\rm s}=\gamma f_{\rm Nyq},$ (2)

and the sampling interval is $T_{\rm s}=1/f_{\rm s}$. Uniform unlimited
sampling of $g(t)$ yields

$\displaystyle y_{n}\triangleq\mathscr{M}_{\lambda}\left(\Re\left\\{g(nT_{\rm
s})\right\\}\right)+{\rm j}\mathscr{M}_{\lambda}\left(\Im\left\\{g(nT_{\rm
s})\right\\}\right),$ (3)

$n=0,1,\cdots,N-1$, where $\mathscr{M}_{\lambda}\left(t\right)$ is the module
mapping defined as

$\displaystyle\mathscr{M}_{\lambda}\left(t\right)\triangleq
2\lambda\left(\left\llbracket\frac{t}{2\lambda}+\frac{1}{2}\right\rrbracket-\frac{1}{2}\right),$
(4)

$\llbracket t\rrbracket\stackrel{{\scriptstyle\text{ def }}}{{=}}t-\lfloor
t\rfloor$ denotes the fractional part of $t$, $\lambda$ denotes the ADC
threshold. Note that the number of samples is $N$. $g_{n}=g(nT_{\rm s})$ can
be written as

$\displaystyle g_{n}=\frac{1}{\sqrt{N}}\sum\limits_{k=1}^{K}c_{k}{\rm e}^{{\rm
j}\omega_{k}n}+w_{n},$ (5)

where $w_{n}\triangleq w(nT_{\rm s})$ denotes the noise sequence, $\omega_{k}$
is

$\displaystyle\omega_{k}=\frac{2\pi}{\gamma}\times\frac{f_{k}}{f_{\rm
max}}\in[0,\frac{1}{\gamma}\times 2\pi].$ (6)

LSE aims to recover $\\{(c_{k},\omega_{k})\\}_{k=1}^{K}$ from
${\mathbf{y}}=[y_{0},\cdots,y_{N-1}]^{\rm{T}}$.

The proposed US-LSE algorithm is motivated by several facts and are presented
below:

  * •

Fact 1 [4]: The original samples can be decomposed as the sum of module
samples and a simple function, i.e., $g_{n}=y_{n}+\epsilon_{g,n}$, and the
first-order finite difference samples also applies.

  * •

Fact 2 [21]: The first-order finite difference of samples of the simple
function are sparse. Besides, their values lie on the grid
$2\lambda{\mathbb{Z}}$.

We make an assumption that the frequencies lie on the grid exactly, i.e.,
$\omega_{k}\in\\{0,\frac{2\pi}{N},\frac{2\pi}{N}\times
2,\cdots,\frac{2\pi}{N}\times(\lfloor\frac{N}{\gamma}\rfloor-1)\\}$. Define
the discrete Fourier transform (DFT) matrix
${\mathbf{F}}\in{\mathbb{C}}^{N\times N}$ with the $(n_{1},n_{2})$ elements
being $\frac{1}{\sqrt{N}}{\rm e}^{-{\rm j}\frac{2\pi(n_{1}-1)(n_{2}-1)}{N}}$,
$n_{1}=1,2,\cdots,N$ and $n_{2}=1,2,\cdots,N$. Obviously,
${\mathbf{F}}{\mathbf{F}}^{\rm H}={\mathbf{F}}^{\rm
H}{\mathbf{F}}={\mathbf{I}}_{N}$. Consequently, $\mathbf{g}$ (5) can be
formulated as

$\displaystyle{\mathbf{g}}={\mathbf{F}}^{\rm H}\bar{\mathbf{c}}+{\mathbf{w}},$
(7)

where $\bar{\mathbf{c}}\in{\mathbb{C}}^{N}$, $\bar{c}_{i}=0$, $i=\lfloor
N/\gamma\rfloor,\lfloor N/\gamma\rfloor+1,\cdots,N$,
$\|\bar{\mathbf{c}}\|_{0}=K$. This implies that the support set of
$\bar{\mathbf{c}}$ is a subset of set $\\{0,1,2,\cdots,\lfloor
N/{\gamma}\rfloor-1\\}$.

Note that

$\displaystyle{\mathbf{g}}={\mathbf{y}}+{\bm{\epsilon}}_{g}.$ (8)

Define the cyclic first-order finite difference
$\Delta{\mathbf{g}}\in{\mathbb{C}}^{N}$ with the $k$th element being

$\displaystyle\Delta g[k]$ $\displaystyle=g_{k+1}-g_{k},k=0,1,\cdots,N-2,$ (9)
$\displaystyle\Delta g[N-1]$ $\displaystyle=g_{0}-g_{N-1},$ (10)

or in matrix form

$\displaystyle\Delta{\mathbf{g}}=({\mathbf{J}}-{\mathbf{I}})\mathbf{g}=({\mathbf{J}}-{\mathbf{I}})({\mathbf{F}}^{\rm
H}\bar{\mathbf{c}}+\mathbf{w}),$ (11)

where ${\mathbf{J}}\in{\mathbb{C}}^{N\times N}$ is a selection matrix which
shifts the row of the identity matrix ${\mathbf{I}}_{N}$ up in one dimension
cyclicly. Applying cyclic first-order finite difference to (8) yields

$\displaystyle\Delta{\mathbf{g}}=\Delta{\mathbf{y}}+\Delta{\bm{\epsilon}}_{g}.$
(12)

Performing DFT on (12) yields

$\displaystyle{\mathbf{F}}\Delta{\mathbf{g}}={\mathbf{F}}\Delta{\mathbf{y}}+{\mathbf{F}}\Delta{\bm{\epsilon}}_{g}.$
(13)

According to
${\mathbf{F}}\Delta{\mathbf{g}}={\mathbf{F}}({\mathbf{J}}-{\mathbf{I}})({\mathbf{F}}^{\rm
H}\bar{\mathbf{c}}+{\mathbf{w}})={\mathbf{d}}\odot(\bar{\mathbf{c}}+\bar{\mathbf{w}})$,
where ${\mathbf{d}}=[1,{\rm e}^{{\rm j}\frac{2\pi}{N}},\cdots,{\rm e}^{{\rm
j}\frac{2\pi(N-1)}{N}}]^{\rm T}-{\mathbf{1}}$, and
$\bar{\mathbf{w}}=\mathbf{F}\mathbf{w}$, one has

$\displaystyle{\mathbf{d}}\odot(\bar{\mathbf{c}}+\bar{\mathbf{w}})=\Delta\tilde{\mathbf{y}}+{\mathbf{F}}\Delta{\bm{\epsilon}}_{g},$
(14)

where $\Delta\tilde{\mathbf{y}}\triangleq{\mathbf{F}}\Delta{\mathbf{y}}$.
According to Fact $2$, $\Delta{\bm{\epsilon}}_{g}$ is a $M$ sparse vector and
$\bar{\mathbf{c}}$ is a $K$ sparse vector. We use least squares to formulate
(14) as the following optimization problem

$\displaystyle\underset{\bar{\mathbf{c}},\Delta{\bm{\epsilon}}_{g}}{\operatorname{minimize}}~{}\|\Delta\tilde{\mathbf{y}}+{\mathbf{F}}\Delta{\bm{\epsilon}}_{g}-\mathbf{d}\odot\bar{\mathbf{c}}\|_{2}^{2},$
(15a) $\displaystyle{\rm subject~{}to}~{}\|\bar{\mathbf{c}}\|_{0}\leq
K,\|\Delta{\bm{\epsilon}}_{g}\|_{0}\leq M,\Delta{\bm{\epsilon}}_{g}\in
2\lambda{\mathbb{Z}}.$ (15b)

Note that $\bar{\mathbf{c}}_{i}=0$, for $i=M+1,M+2,M+3,\cdots,N$. Denote
$\Delta\tilde{\mathbf{y}}_{\rm sub}$ as a subvector of
$\Delta\tilde{\mathbf{y}}$ with index chosen from $M+1$ to $N$,
${\mathbf{F}}_{\rm sub}$ denotes a submatrix chosen from the rows of
${\mathbf{F}}$ with index chosen from $M+1$ to $N$. Then we solve the
subproblem

$\displaystyle\underset{\Delta{\bm{\epsilon}}_{g}}{\operatorname{minimize}}~{}\|\Delta\tilde{\mathbf{y}}_{\rm
sub}+{\mathbf{F}}_{\rm sub}\Delta{\bm{\epsilon}}_{g}\|_{2}^{2},$ (16a)
$\displaystyle{\rm subject~{}to}~{}\|\Delta{\bm{\epsilon}}_{g}\|_{0}\leq P.$
(16b)

via OMP. A sufficient condition for recovering $\Delta{\bm{\epsilon}}_{g}$ is
that ${\rm spark}({\mathbf{F}}_{\rm sub})>2P$ [23, Corollary 1], where ${\rm
spark}({\mathbf{A}})$ denotes the smallest number of columns of $\mathbf{A}$
that are linearly dependent. Obviously, ${\rm spark}({\mathbf{F}}_{\rm
sub})\in[2,N-M]$. Therefore, $N-M=N(1-\frac{1}{\gamma})\geq 2P$ should be
satisfied. This implies that the oversampling factor should satisfy

$\displaystyle{\gamma}\geq\frac{1}{1-2P/N}.$ (17)

Once $\Delta\bm{\epsilon}_{g}$ is recovered, $\bar{\mathbf{c}}$ can be
estimated via (15). Overall the US-LSE is summarized in Algorithm 1.

Algorithm 1 LSE via Unlimited Sampling (US-LSE)

1: Input: Modulo samples ${\mathbf{y}}\in{\mathbb{C}}^{N}$, number of
sinusoids $K$, and the number of folding instants $P$.

2: Compute $\Delta{\mathbf{y}}$ and
$\Delta\tilde{\mathbf{y}}\triangleq{\mathbf{F}}\Delta{\mathbf{y}}$.

3: Solve (16) to obtain $\widehat{\Delta{\bm{\epsilon}}}_{g}$.

4: Round $\widehat{\Delta{\bm{\epsilon}}}_{g}$ to the nearest grid
$2\lambda{\mathbb{Z}}$.

5: Solve (15) and return $\widehat{\bar{\mathbf{c}}}$.

##  3 Numerical Simulation

In this section, numerical experiments are conducted to verify the performance
of the proposed US-LSE. The parameters are set as follows: $N=256$, the number
of frequencies is $K=3$, the frequencies are randomly generated from
$\\{\frac{2\pi}{N},2\frac{2\pi}{N}\cdots,\lfloor\frac{N}{\gamma}\rfloor\frac{2\pi}{N}\\}$,
$\lambda=0.06$, the amplitudes are identical and phases are randomly drawed
from $[0,2\pi)$. The SNR is defined as
$\rm{SNR}=\|\mathbf{x}\|^{2}_{2}/\mathbb{E}[\|\mathbf{w}\|^{2}_{2}]$, where
$\mathbf{x}$ is the noise-free signals and $\mathbf{w}$ is the noise. The
normalized mean square error (NMSE) is
$\rm{NMSE}(\widehat{\mathbf{x}})=\|\widehat{\mathbf{x}}-{\mathbf{x}}\|^{2}_{2}/\|{\mathbf{x}}\|^{2}_{2}$,
where $\widehat{\mathbf{x}}$ denotes the estimate of the algorithm. The US-LSE
is compared with the existing method named USE that recovers the spectrum from
modulo-folded samples [5, 4]. We define the event in which the algorithm
successfully recovers the true signal $\mathbf{x}$ provided that
$\rm{NMSE}(\widehat{\mathbf{x}})<-15$ dB. All the results are averaged over
$300$ Monte Carlo (MC) trials.

The probabilities of successful recovery of USE and US-LSE algorithms versus
the SNR and oversampling factor $\gamma$ are shown in Fig. 1. With the
oversampling factor $\gamma$ being fixed and $\gamma\leq 10$, when
$\rm{SNR}\leq 25$ dB, the probability of successful recovery of USE is lower
than that of US-LSE, when $\rm{SNR}\geq 30$ dB, the probability of successful
recovery of USE is higher than that of US-LSE. This demonstrates that in the
low oversampling scenario where $\gamma\leq 10$, US-LSE is more stable in the
lower SNR scenario where $\rm{SNR}\leq 25$ dB, while USE performs better than
US-LSE in the higher SNR scenario where $\rm{SNR}\geq 30$ dB. The poor
performance of USE is due to the high order difference operator which may
amplify noise in the low SNR scenario. For $\gamma\geq 12$, the phenomenon is
similar to the scenario in which $\gamma\leq 10$ when the SNR is low, while
both algorithms successfully recover the true signal in the high SNR scenario.
This demonstrates that US-LSE benefits more from oversampling, as oversampling
increases the number of effective observations $\Delta\tilde{\mathbf{y}}_{\rm
sub}$ and $\Delta{\bm{\epsilon}}_{g}$ can be recovered more easily.

Fig. 1: The probability of successful recovery versus SNR and oversampling
factor $\gamma$.

##  4 Empirical Experiments

Empirical experiments are conducted to demonstrate the performance of US-LSE
by using an AWR1642 radar. The sampling frequency is $f_{s}=5$ MHz. The chirp
rate is $\kappa=15\times 10^{12}$ Hz/s. The maximum distance is
$r_{\rm{max}}=cf_{s}/(2\kappa)=50$ m. Provided that the distances of targets
are less than $5$ m which is the case in our experiments, it can be shown that
the range mesearment model can be formulated as an oversampling line spectra
estimation model with $\gamma=10$. The USF measurements are generated in
software from the original measurements.

The number of measurements is $N=256$, the number of targets is $K=2$ and
$M=25$. The maximum amplitude of the inphase and quadrature components is
about $240$, the threshold is $\lambda=80$. We conduct two experiments, as
shown in Fig. 2. We run the Newtonalized OMP (NOMP) to process the original
measurements and output the range estimates as the ground truth [3]. Note that
US-LSE is the on-grid approach and the grid spacing is $r_{\rm{max}}/N\approx
0.195$m. For scenario 1, the radial distances of two corners named corner 1
and corner 2 and $1.421$m and $2.171$m, respectively. The grids closest to
corner 1 and corner 2 are $1.367$m and $2.148$m, and the off-grid distances
are $0.054$m and $0.023$m. For scenario 2, the radial distances of corner 1
and corner 2 and $1.421$m and $2.419$m, respectively. The grids closest to
corner 1 and corner 2 are $1.367$m and $2.344$m, and the off-grid distances
are $0.054$m and $0.075$m. We estimate the SNRs of the two targets based on
the spectrum of the original data, which are $11.5$ dB and $24.3$ dB for
corner 1 and corner 2 in scenario 1, $14.0$ dB and $22.0$ dB for corner 1 and
corner 2 in scenario 2.

For scenario 1, we apply US-LSE to perform LSE. Note that the number of folds
of the modulo non-linearity is $P=85$ satisfying (17). Fig. 3 shows the
recovery results of US-LSE and USE. USE failed to recover the inpahse
component and recover the quadrature component, while US-LSE successfully
recovers both the inphase and quadrature components. For scenario 2, Fig. 3
shows that USE failed to recover both the inphase and quadrature components,
while US-LSE still recovers the signal.

Once the $\Delta\bm{\epsilon}_{g}$ is estimated, we solve (15) to obtain
$\tilde{\bar{\mathbf{c}}}$. The results are shown in Fig. 5. For scenario 1,
the recovery results of US-LSE match the ground truth quite well. For scenario
2, the amplitude of corner 1 is consistent with ground truth, and the
amplitudes of two bins closest to the distance of target 2 are large, due to
the significant off-grid effects of corner 2.

Fig. 2: Field experiments setup. (a) The radial distances of corner $1$ and
corner $2$ are $1.42$m and $2.17$m, respectively. (b) The radial distance of
corner $1$ and corner $2$ are $1.42$m and $2.42$m. Fig. 3: The recovery
results of scenario $1$. (a) The inphase component, (b) the recovered inphase
component, (c) the quadrature component, (b) the recovered quadrature
component. Fig. 4: The recovery results of scenario $2$. (a) The inphase
component, (b) the recovered inphase component, (c) the quadrature component,
(b) the recovered quadrature component. Fig. 5: The estimated
$\tilde{\bar{\mathbf{c}}}$ for scenario $1$ and scenario $2$.

##  5 Conclusion

In this paper, the LSE is studied in the context of USF. The US-LSE is
proposed, which exploits the oversampling and the bandlimited property of the
line spectrum to recover the folding instants and the line spectral. The
robustness of US-LSE compared to the existing approach is demonstrated through
numerical experiments. Besides, we also evaluate the performance of US-LSE by
processing the data acquired by mmWave radar. It is shown that the US-LSE
successfully recovers the LSE with the oversampling factor $\gamma=10$.
Several relevant future concerns are worth noting. For example, we recover the
first-order difference of the simple functions by solving a relaxed problem.
It may further improve the robustness and performance of the proposed approach
if one could jointly estimate the frequencies and the first-order difference
of the simple functions.

## References

  * [1] P. Stoica, R. L. Moses, et al., Spectral analysis of signals, vol. 452.  Pearson Prentice Hall Upper Saddle River, NJ, 2005.
  * [2] T. L. Hansen, B. H. Fleury, and B. D. Rao, “Superfast line spectral estimation,” IEEE Transactions on Signal Processing, vol. 66, no. 10, pp. 2511–2526, 2018.
  * [3] B. Mamandipoor, D. Ramasamy, and U. Madhow, “Newtonized orthogonal matching pursuit: Frequency estimation over the continuum,” IEEE Transactions on Signal Processing, vol. 64, no. 19, pp. 5066–5081, 2016.
  * [4] A. Bhandari, F. Krahmer, and R. Raskar, “On unlimited sampling and reconstruction,” IEEE Transactions on Signal Processing, vol. 69, pp. 3827–3839, 2020.
  * [5] A. Bhandari, F. Krahmer, and R. Raskar, “Unlimited sampling of sparse sinusoidal mixtures,” in 2018 IEEE International Symposium on Information Theory (ISIT), pp. 336–340, IEEE, 2018.
  * [6] S. Fernández-Menduina, F. Krahmer, G. Leus, and A. Bhandari, “Doa estimation via unlimited sensing,” in 2020 28th European Signal Processing Conference (EUSIPCO), pp. 1866–1870, IEEE, 2021.
  * [7] S. Fernández-Menduiña, F. Krahmer, G. Leus, and A. Bhandari, “Computational array signal processing via modulo non-linearities,” IEEE Transactions on Signal Processing, vol. 70, pp. 2168–2179, 2022.
  * [8] S. Rudresh, A. Adiga, B. A. Shenoy, and C. S. Seelamantula, “Wavelet-based reconstruction for unlimited sampling,” in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 4584–4588, IEEE, 2018.
  * [9] T. Feuillen, M. Alaee-Kerahroodi, A. Bhandari, B. Ottersten, et al., “Unlimited sampling for fmcw radars: A proof of concept,” in 2022 IEEE Radar Conference (RadarConf22), pp. 1–5, IEEE, 2022.
  * [10] V. Bouis, F. Krahmer, and A. Bhandari, “Multidimensional unlimited sampling: A geometrical perspective,” in 2020 28th European Signal Processing Conference (EUSIPCO), pp. 2314–2318, IEEE, 2021.
  * [11] O. Graf, A. Bhandari, and F. Krahmer, “One-bit unlimited sampling,” in ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5102–5106, IEEE, 2019.
  * [12] A. Bhandari, O. Graf, F. Krahmer, and A. I. Zayed, “One-bit sampling in fractional fourier domain,” in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 9140–9144, IEEE, 2020.
  * [13] N. I. Bernardo, J. Zhu, Y. C. Eldar, and J. Evans, “Design and analysis of hardware-limited non-uniform task-based quantizers,” arXiv preprint arXiv:2208.07525, 2022.
  * [14] A. Bhandari, “Back in the us-sr: Unlimited sampling and sparse super-resolution with its hardware validation,” IEEE Signal Processing Letters, vol. 29, pp. 1047–1051, 2022.
  * [15] A. Bhandari and F. Krahmer, “Hdr imaging from quantization noise,” in 2020 IEEE International Conference on Image Processing (ICIP), pp. 101–105, IEEE, 2020.
  * [16] D. Florescu, F. Krahmer, and A. Bhandari, “The surprising benefits of hysteresis in unlimited sampling: Theory, algorithms and experiments,” IEEE Transactions on Signal Processing, vol. 70, pp. 616–630, 2022.
  * [17] D. Florescu and A. Bhandari, “Unlimited sampling with local averages,” in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5742–5746, IEEE, 2022.
  * [18] D. Florescu and A. Bhandari, “Unlimited sampling via generalized thresholding,” in 2022 IEEE International Symposium on Information Theory (ISIT), pp. 1606–1611, IEEE, 2022.
  * [19] O. Musa, P. Jung, and N. Goertz, “Generalized approximate message passing for unlimited sampling of sparse signals,” in 2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pp. 336–340, IEEE, 2018.
  * [20] G. Lu and H. Liu, “High dynamic range sensing using multi-channel modulo samplers,” in 2020 IEEE 11th Sensor Array and Multichannel Signal Processing Workshop (SAM), pp. 1–5, IEEE, 2020.
  * [21] A. Bhandari, F. Krahmer, and T. Poskitt, “Unlimited sampling from theory to practice: Fourier-prony recovery and prototype adc,” IEEE Transactions on Signal Processing, vol. 70, pp. 1131–1141, 2021.
  * [22] E. Azar, S. Mulleti, and Y. C. Eldar, “Residual recovery algorithm for modulo sampling,” in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 5722–5726, IEEE, 2022\.
  * [23] D. L. Donoho and M. Elad, “Optimally sparse representation in general nonorthogonal dictionaries via l1 minimization,” Proceedings of the National Academy of Sciences, vol. 100, no. 5, pp. 2197–2202, 2003.