---
abstract: 'Quantum criticality beyond the Landau paradigm represents a fundamental problem in condensed matter and statistical physics. Heavy fermion systems with multipolar degrees of freedom can play an important role in the search for its universal description. We consider a Kondo lattice model with both spin and quadrupole degrees of freedom, which we show to exhibit an antiferroquadrupolar phase. Using a field theoretical representation of the model, we find that Kondo couplings are exactly marginal in the renormalization group sense in this phase. This contrasts with the relevant nature of the Kondo couplings in the paramagnetic phase and, as such, it implies that a Kondo destruction and a concomitant small to large Fermi surface jump must occur as the system is tuned from the antiferroquadrupolar ordered to the paramagnetic phase. Implications of our results for multipolar heavy fermion physics in particular and metallic quantum criticality in general are discussed.'
author:
- 'Hsin-Hua Lai'
- 'Emilian M. Nica'
- 'Wen-Jun Hu'
- 'Shou-Shu Gong'
- Silke Paschen
- Qimiao Si
bibliography:
- 'biblio4AFQ.bib'
title: 'Kondo Destruction and Multipolar Order – Implications for Heavy Fermion Quantum Criticality'
---

*Introduction*— In strongly correlated systems, multiple building blocks often interplay with each other and create a variety of quantum phases and their transitions. Examples include the spin, orbital and nematic degrees of freedom in the iron-based systems [@Kamihara2008; @NatRevMat:2016], which lead to a rich landscape of electronic orders, and the spin and valley degrees of freedom in bilayer graphenes twisted by magic angles [@YCao1; @YCao2], which appear to yield a surprising Mott insulator near which superconductivity develops. The multiple degrees of freedom allow for not only the commonly observed antiferromagnetic (AF) states, but also “hidden" orders, with unusual order parameters that cannot readily be probed by experiments directly. A prominent example is the quadrupolar order, which breaks the spin-rotational symmetry as in any conventional magnetic order but, unlike the latter, preserves the time-reversal symmetry. Such an order has been proposed for frustrated magnetic systems [@Papanicolaou1988; @Tsunetsugu2006; @Lauchli2006; @Smerald2013] and even for the nematic phase of the iron-chalcogenide FeSe [@YuSi_AFQ; @QisiWang2016; @Lai_AFQ2017]. Multipolar degrees of freedom are also being discussed in noncollinear antiferromangets [@Suzuki:2017]. They also arise in many heavy fermion metals, producing a variety of fascinating properties [@Custers2012; @Martelli2017; @Sakai2011; @Lee2018; @Bauer2002; @McCollam2013].

Heavy fermion compounds typically involve local spin moments, which experience RKKY interactions between each other and Kondo interactions with conduction electrons, and exhibit quantum phase transitions between paramagnetic and AF ground states [@SiSteglich2010; @ColemanSchofield2005]. While the Kondo effect has been a hallmark of heavy fermion physics, a Kondo destruction has been shown to arise from the dynamical competition between the RKKY and Kondo interactions  [@Si-Nature; @Colemanetal]. It has been demonstrated in studies of Kondo lattice models from both the paramagnetic [@Si-Nature] and AF-ordered  [@Yamamoto2007] sides. Because the Kondo destruction yields quantum criticality that is beyond the Landau framework of order-parameter fluctuations, it is important to assess its universality by considering settings that involve other types of local degrees of freedom.

An especially opportune setting arises in heavy fermion systems with co-existing local spin and multipolar moments [@Custers2012; @Martelli2017], which allow for not only AF orders but also quadrupolar ones. In Ce$_3$Pd$_{20}$Si$_6$, an antiferroquadrupolar (AFQ) order has been experimentally determined [@Por16.1], and a sequence of quantum critical points was discovered upon tuning by a magnetic field [@Martelli2017]. Theoretical calculations that approach the transitions from the paramagnetic side demonstrated a sequential Kondo destruction [@Martelli2017]. This provides the motivation to study the Kondo effect and Fermi surface in the AFQ ordered state.

In this Letter, we address this pressing problem using a spin-$1$ Kondo lattice model, which contains both spin and quadrupole local degrees of freedom. We demonstrate a robust AFQ phase, and describes its low-energy effective theory in terms of a quantum non-linear sigma model (NLsM) [@Smerald2013]. Adapting a combined boson-fermion renormalization group (RG) procedure [@Yamamoto2010], we show that the Kondo couplings are exactly marginal in the RG sense and thereby establish a Kondo destruction in the multipolar order.

The model we consider is $H_{KL} = H_S + H_c + H_K$, with $$\begin{aligned}
&& H_{S} = \sum_{i j}\left[ J_{ij}\left( \bm{S}_i \cdot \bm{S}_j \right) + K_{ij} \left(\bm{S}_i \cdot \bm{S}_j\right)^2\right], \label{Eq:BBQ}\\ 
&& H_c = \sum_{\bm{k}, \alpha = x,y,z} \epsilon_{\bm{k}} \psi_{\bm{k} \alpha}^\dagger \psi_{\bm{k} \alpha}, \\
 && H_K = \sum_j \left( J_K^I \bm{S}_j \cdot \bm{s}_{c,j} + J_K^{II} \bm{Q}_j \cdot \bm{q}_{c,j}\right), \label{Eq:Kondo}\end{aligned}$$ where $H_S$ represents the spin-$1$ bilinear-biquadratic Hamiltonian, in which we choose $J_{ij} = J_n,~K_{ij}= K_n $ for $i, j$ connected by $n$th-neighbor bonds. The spin-$1$ nature implies that the existence of local quadrupolar moments. The 5-component quadrupolar operator at site $i$, ${\bf Q}_i$, can be defined as: $Q^{x^2 - y^2}_i = (S^x_i)^2 - (S^y_i)^2$, $Q^{3z^2 - r^2}_i = [ 2 (S^z_i)^2 - (S^x_i)^2 - (S^y_i)^2]/\sqrt{3}$, $Q^{xy}_i = S^x_i S^y_i + S^y_i S^x_i$, $Q^{yz}_i =  S^y_iS^z_i +S^z_i S^y_i$, and $Q^{zx}_i = S^z_i S^x_i + S^x_i S^z_i$. The biquadratic term can be re-expressed as $({\bf S}_i \cdot {\bf S}_j )^2 = ({\bf Q}_i \cdot {\bf Q}_j) /2 
- ({\bf S}_i \cdot {\bf S}_j) /2 + ({\bf S}^2_i {\bf S}^2_j)/3$. At the high-symmetry point, $J_n = K_n$, the symmetry is enhanced from SU(2) to SU(3). Here, the spin and quadrupolar moments can be transformed to each other under SU(3) rotations. Our focus will be on the AFQ phase away from the SU(3) point; however, as we will see, the time-reversal-invariant basis that is natural for the SU(3) point – which can be related to the $s^z = \pm 1, 0$ basis under a unitary transformation – will greatly facilitate our analysis. $H_c$ describes the conduction electrons, which have three flavors with flavor index $\alpha = x, y, z$ in the SU(3) time-reversal-invariant basis. Within the 3-flavor conduction electron description, both the electrons’ spin $\bm{s}_c$ and their $5$-component quadrupoles, $\bm{q}_c$, are expressed in bilinear forms (see Supplemental Materials for explicit forms of $\bm{q}_c$). $H_K$ represents the Kondo couplings between the local moments and conduction electrons.

*Existence of $(\pi,\pi)$ antiferroquadrupolar order*— We first study the spin-$1$ bilinear-biquadratic lattice model, Eq. , numerically using the large-scale Density Matrix Renormalization Group (DMRG) analysis. It is known that this SU(3)-symmetric point can host a phase with both spin and quadrupolar orders, which can be transformed to each other under SU(3) transformations [@Papanicolaou1988; @Tsunetsugu2006; @Lauchli2006; @Smerald2013]. Away from the SU(3) point, we find that increasing the weight of the biquadratic terms can stabilize the quadrupolar order. To illustrate the robustness of the quadrupolar order in this model, we fix $J_1=1$, $K_1 = 1.2$, and determine the phase diagram as a function of $K_2/J_1 (J_2 = 0)$, which is shown in Fig. \[Fig:DMRG\_result\]. We find the the AFQ order with $\bm{q}=(\pi,\pi)$. This $(\pi,\pi)$-AFQ is a two-sublattice order characterized by the staggered expectation values of the $Q^{x^2 - y^2}$, *i.e.*, $\langle Q^{x^2 - y^2}_j \rangle \sim (-1)^{j}$. The $(\pi,\pi)$-AFQ phase can be identified by calculating the spin structure factor ($m^2_{\bm{S}}$) and quadrupolar structure factor ($m^2_{\bm{Q}}$). As illustration, we consider the case $K_2/J_1 = - 0.3$, which shows strong peak at $\bm{q} = (\pi,\pi)$ for $m^2_{\bm{Q}}$ and much weaker peak at $\bm{q} = (\pi,0)/(0,\pi)$ for $m^2_{\bm{S}}$, shown in the bottom left panel of Fig. \[Fig:DMRG\_result\] for system size $L_y = 8$. Performing the finite-size scaling analysis, bottom right panel of Fig. \[Fig:DMRG\_result\], we find nonzero $m^2_{\bm{Q}}$ and vanishing $m^2_{\bm{S}}$, which shows the presence of the $(\pi,\pi)$-AFQ order.

         

*Non-linear sigma model for $(\pi,\pi)$-AFQ*— Because the commensurate AFQ breaks the spin-rotational symmetry but is time-reversal invariant, we can expect three Goldstone modes. To specify the low-energy effective theory including couplings that involve the conduction electrons, we describe the Kondo lattice model using a NLsM representation by adapting the method illustrated in Ref. [@Smerald2013]. We first introduce the SU(3) time-reversal-invariant basis, $ |x\ra = i ( | 1 \ra - |\bar{1})/\sqrt{2},~~~ |y\ra = ( | 1\ra + | \bar{1} )/\sqrt{2}, ~~~|z \ra = -i | 0 \ra,$ where $|S^z = 1\ra \equiv |1\ra$. The state at site $j$ can be written in terms of $ {\bf d}_j = (d^x_j,~d^y_j,~d^z_j)$, where $d_j^\alpha$ are complex numbers, with constraints from the normalization and from fixing the global phase among $d^{\alpha=x,y,a}$, *i.e.*, ${\bf d}_j \cdot \bar{\bf d}_j = 1$, and ${\bf d}_j^2 = \bar{\bf d}_j^2$, where $\bar{\bm{d}}$ means complex conjugate of ${\bm{d}}$. The Hamiltonian is then re-expressed as (see Supplemental Materials [@supp]) $$\begin{aligned}
H_{S} &=& \sum_{i j } \left[ J_n |{\bf d}_i \cdot \bar{\bf d}_j|^2 + J_n' | {\bf d}_i \cdot {\bf d}_j |^2  \right] \equiv H_{SU(3)} + H',~~~\end{aligned}$$ where we have ignored the inconsequential constant terms. We define the deviations from the SU(3) couplings, $J_n' \equiv K_n - J_n$, and explicitly separate the SU(3)-invariant part of the Hamiltonian $H_{SU(3)}$ from the SU(3)-breaking part $H'$. At the SU(3) point, *i.e.,* $J_n'=0$ or $J_n = K_n$, we can see that $H_S$ becomes a pure function of $|(\bm{d}_i \cdot \bar{\bm{d}}_j)|^2$. For the lowest-energy real-space pattern, we need to minimize the nearest-neighbor $|\bm{d}_i\cdot \bar{\bm{d}}_j|$ and maximize the $2$nd-neighbor $|\bm{d}_i \cdot \bar{\bm{d}}_j |$. For obtaining the NLsM description, we choose the ground state pattern of the $(\pi,\pi)$-AFQ, which satisfies the above requirement, as $ ({\bf d}^{gs}_A)^T = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}$, and $({\bf d}^{gs}_B)^T = \begin{pmatrix} 0 & 1 & 0 \end{pmatrix}$. Such a pattern is illustrated in Fig. \[Fig:pipiAFQ\], which gives the correct $(\pi,\pi)$-AFQ order at the semi-classical level.

Starting from the SU(3) point, we know that the ground state energy is invariant under the global rotation ${\bf d} \rightarrow {\bf U}{\bf d}$, provided that ${\bf U}^{-1} = {\bf U}^\dagger$. To describe global rotations, we find that Gell-Mann matrices provide a natural choice of basis at the SU(3)-symmetric point. In general, we require $6$ distinct generators for SU(3). However, for the present $(\pi,\pi)$-AFQ phase, only $4$ out of $6$ are needed. The $4$ Gell-Mann matrices we choose are represented by $\lambda_{1- 4}$ (see Supplemental Materials for the explicit matrix forms). The global rotations in the complex space can be expressed as ${\bf U}(\phi) \equiv exp[{i\sum_{p =1}^4 \lambda_p \phi_p}].$ Besides the global rotations that preserve the ground state energy, we also need to consider the rotations involving the canting of the directors of ground state configurations, which increase the energy. The canting fields are represented as $\mu_{1 \sim 4}$ (see Supplemental Materials for details). The general rotations can be represented as $$\begin{aligned}
\label{Eq:rotation_m}
\nonumber {\bf D}(\phi, \ell, v) & = &  e^{i \sum_{p=1}^4 \lambda_p \phi_p + i \mu_1 \ell^z_1 + i \mu_2 v_A + i \mu_3 v_B + i \mu_4 \ell_2^z}\\
 && \hspace{-1.5 cm} \simeq U(\phi) \sum_{n=0}^\infty \left( i \mu_1 \ell^z_1 + i \mu_2 v_A + i \mu_3 v_B + i \mu_4 \ell_2^z\right)^n ,\end{aligned}$$ where we approximately separate out the global rotation matrix $U(\phi)$ and Taylor expand the terms involving the canting fields. The general configuration of ${\bf d}_{a=A/B}$ can be obtained by applying the general rotation on the chosen ground state configurations ${\bf d}^{gs}_{a=A/B}$, *i.e.*, $\bm{d}_a = \bm{D}(\phi,\ell,v) \bm{d}^{gs}_a$. To obtain the low-energy descriptions within the harmonic theory, we keep the Taylor expansion up to $n=2$. Introducing $\ell^z \equiv \ell_1^z + i \ell^z_2$ and re-parametrizing $$\begin{aligned}
\label{Eq:U(n)}
{\bf U} = \begin{pmatrix}
n^x_A & n^x_B & n^x_C \\
n^y_A & n^y_B & n^y_C \\
n^z_A & n^z_B & n^z_C 
\end{pmatrix} = \begin{pmatrix}
\bm{n}_A & \bm{n}_B & \bm{n}_C
\end{pmatrix},\end{aligned}$$ where the vectors ${\bf n}^T_a = \begin{pmatrix} n_a^x & n_a^y & n_a^z \end{pmatrix}$ inherit the constraints of ${\bf d}$ with ${\bf n}_a \cdot \bar{\bf n}_{b\not=a} = 0$, ${\bf n}_a^2 = {\bf n}_a^2$, and ${\bf n}_a \cdot \bar{\bf n}_a = 1$, with the vector ${\bf n}_c \equiv \bar{\bf n}_A \times \bar{\bf n}_B$ being introduced as a convenient piece of book-keeping in the present $(\pi, \pi)$-AFQ. It is not an independent degrees of freedom, we can fully re-express $\bm{d}_{a = A/B}$ as functions of $\bm{n}_a$, $\ell^z$, and $v_{A/B}$ (see Supplemental Materials). Taylor expanding $U(\phi)$ and keeping only the leading linear terms in $\phi$, we can see that $\bm{n}_a$ and the bosonic field $\phi_{1\sim 4}$ are related by $$\begin{aligned}
 \label{Eq:n_to_phi}
\bm{n}_A \simeq \begin{pmatrix}
1 \\
-\phi_1 + i \phi_4 \\
\phi_2
\end{pmatrix};&~~
\bm{n}_B \simeq \begin{pmatrix}
\phi_1 + i \phi_4 \\
1 \\
-\phi_3
\end{pmatrix},\end{aligned}$$ where $\bm{n}_C = \bar{\bm{n}}_A \times \bar{\bm{n}}_B$. Away from the SU(3) point, we assume that $H'$ can be treated perturbatively and does not affect our results in any significant manner.

Adopting the strategy of Ref. [@Smerald2013], we partition the square lattice into clusters (Figs. \[Fig:cluster\_network\]-\[Fig:cluster\_single\]), each of which containing 8 bonds \[4 nearest-neighbor bonds and 4 second-neighbor bonds (dashed lines)\], and perform the gradient expansion about the center of a cluster. Within the cluster picture, the partition function can be concisely expressed as $\mathcal{Z}_\square = \int D[{\bf d}] e^{- S_\square [{\bf d}]}$, where the action includes the kinetic terms $S_{kin}$ and the Hamiltonian terms $S_H$, $S_\square = S_{kin} + S_{Hs}$. The continuous descriptions of $S_{kin}=\int_0^\beta d\tau \mathcal{L}_{kin}$ and $S_{Hs}=\int_0^\beta d\tau H_S = \int_0^\beta d\tau \left[ H_{SU(3)} + H' \right]$ can be straightforwardly obtained after the gradient expansion. The detailed results are presented in the Supplemental Materials, and all the terms ($\mathcal{L}_{kin},~H_s$) are functions of $\ell_1^z$, $\ell^z_2$, $v_{A/B}$, and $\bm{n}_{a} $. Integrating out the canting fields by solving the differential equations, $
\delta \mathcal{L}/\delta \ell_1^z = 0$, $\delta \mathcal{L}/ \delta \ell_2^z =0$, $\delta \mathcal{L}/\delta v_A = 0$, $\delta \mathcal{L}/ \delta v_B =0$ within the steepest-descent approximation, we obtain the NLsM for the $(\pi,\pi)$-AFQ at the harmonic level. Stability requires $K_1>0$, $ K_2 <0$, $K_1 - J_1 \geq 0$, $J_2 - K_2 >0$, $J_1 + J_2 - K_2 > 0$, and $K_1 + 2K_2 - 4J_2 \geq 0$. Focusing on the regime away from the SU(3) point, $J_n' \not=0$, we find that $\phi_4$ can be ignored since it represents the spin-wave mode that is always gapped due to the finite mass term. The effective NLsM for the $(\pi,\pi)$-AFQ is $$\begin{aligned}
\label{Eq:Ls}
\nonumber \mathcal{L}^{NLsM}_S \simeq && \frac{ ( \partial_\tau \phi_1 )^2}{4(J_1 + J_2 - K_2)} + 2 (K_1 - 2 K_2) \sum_\lambda (\partial_\lambda \phi_1)^2  \\
&& + \sum_{a=2,3} \left[ \frac{(\partial_\tau \phi_a)^2}{8( J_2 - K_2)} - 2K_2 \sum_\lambda (\partial_\lambda \phi_a)^2\right]. \end{aligned}$$ *Kondo couplings*— Using the identity $\bm{d}_a = D(\phi,\ell,v) \bm{d}^{gs}_a$, and Eqs. -, we can straightforwardly write down the fluctuating $8$-component spin/quadrupolar field in the NLsM description (detailed in Supplemental Materials). Concisely, the $8$-component field can be separate into a uniform part and an oscillating part, $$\begin{aligned}
Q(\bm{r}) \simeq \mathcal{Q}_{\bm{0}} + \mathcal{Q}_{\bm{M}} \cos(\bm{M}\cdot \bm{r}),~~~~~\end{aligned}$$ where $\bm{M} \equiv (\pi,\pi)$, $\mathcal{Q}_{\bm{q}}$ represents the low-energy $8$-component field with momenta $\bm{q}=\bm{0}, \bm{M}$. We remark that that the uniform part contains a “static" background of $Q^{3z^2 -r^2}$ that can directly couple to the $q_c^{3z^2 -r^2}$ of the $3$-flavor conduction electrons due to he Kondo couplings, Eq. . This static background $Q^{3z^2 - r^2}$ field is only invariant under rotation between $x$-$y$ plane, which breaks the SU(3) symmetry of the conduction elections down to SU(2)$\times$U(1), where SU(2) is spanned by the $c_x$ and $c_y$ conduction elections and U(1) is spanned by $c_z$. We expect that the Fermi velocities of $c_x$ and $c_y$ are the same ($v_x = v_y$) but is different from that of $c_z~(v_z)$.

To be specific, at low energies, the conduction electrons in the presence of the static quadrupolar background is $$\begin{aligned}
\label{Eq:Sc}
\nonumber S_c && \simeq  \int d^d\bm{K} d\epsilon \sum_{\alpha = x,y} \psi^\dagger_\alpha (\bm{K}, i \epsilon) ( i \epsilon - \xi_K)\psi_\alpha(\bm{K},i\epsilon) + \\
&& + \int d^d\bm{K}' d\epsilon' \psi^\dagger_z(\bm{K}',i\epsilon')(i\epsilon' - \tilde{\xi}_{K'}) \psi_z(\bm{K}',i\epsilon'),\end{aligned}$$ where $\xi_K = v_F(K - K_F)$, and $\tilde{\xi}_{K'} = \tilde{v}_{F}(K' - K'_F)$, where $K_F$ and $K'_F$ are Fermi momenta for $c_{x/y}$ and $c_z$ and are generically different. The spin dipolar and quadrupolar degrees of freedom consisting of $\psi_{x/y}$ and $\psi_{z}$ fermions, $e.g.$, $s^x_c \sim \psi^\dagger_y \psi_z - \psi^\dagger_z \psi_y$, can be ignored due to the finite energy gap $(\Delta E)$ between the $\psi_{x/y}$ bands and the $\psi_{z}$ band, $\Delta E \propto J_K \la Q^{3z^2 - r^2} \ra$.

Shifting our focus to the Kondo couplings, Eq. , we can now re-express them as $$\begin{aligned}
\nonumber \mathcal{L}_K = && -J_K^I   \left( 2\ell^z_2 + \frac{i}{2} \left(n^y_A - \bar{n}^y_A - n^x_B  + \bar{n}^x_B\right)\right)s_c^z - \\
&&  - J_K^{II}  \left( 2 \ell^z_1 + \frac{1}{2}\left( n^y_A + \bar{n}^y_A + n^x_B + \bar{n}^x_B \right)\right) q_c^{xy},~~~\end{aligned}$$ where only the uniform part $\mathcal{Q}_{\bm{0}}$ couples to the conduction electrons near the Fermi surface. We integrate out the canting field within the steepest descent approximation. After some algebra, we conclude that the effective low-energy description of the action is $S = S^{NLsM}_s + S_K + S_c$, where $S^{NLsM}_S$ and $S_c$ are defined in Eqs.  and , and $S_K$ is $$\begin{aligned}
S_K & \simeq & \int_0^\beta \frac{d\tau}{2} \int d^2 r  \lambda_z \left( s^z_c \partial_\tau \phi_1 \right),\end{aligned}$$ where $\lambda_z = -i J_K^I/[4(J_1 + J_2 - K_2 )]$ is the dimensionless coupling.

*Exact marginality of Kondo couplings*— We now analyze the scaling of the Kondo coupling $\lambda_z$ in the $(\pi,\pi)$-AFQ using the RG procedure described in Ref. [@Yamamoto2010]. For clarity, we introduce $\Phi_a(\bm{r},\tau) \equiv \partial_\tau \phi_a(\bm{r},\tau)$. The scaling dimension of $\Phi_a(\bm{r},\tau)$ can be directly read out, $\Delta[\Phi_a(\bm{r},\tau)] = 1$, indicating that the scaling dimension of its Fourier partner as $\Delta[\Phi_a (\bm{q},\omega)] = -d$, where $d$ is the spatial dimension. For the conduction electron fields, we obtain that $\Delta[\psi_c(\bm{K}, \omega)] = -3/2$. We can see that, at the tree level, the Kondo coupling is *marginal*, $\Delta[S_K] = \Delta[ dk d\epsilon d^dq d \omega \psi^{x \dagger}_c (\bm{k} + \bm{q}, \epsilon + \omega) \psi^y_c (\bm{k}, \omega) \Phi_1 (\bm{q},\omega)] = 1 + 1 + d + 1 + 2 (-3/2) -d = 0$ (see Supplemental Materials). .

We then turn to what happens beyond the tree level. Considering a spherical Fermi surface of conduction electrons, we approximate their contribution via a momentum integral near Fermi surface. Keeping the most relevant term, we obtain $\int d^d \bm{K} = \int^{K_F + \Lambda}_{K_F - \Lambda} K^{d-1} dK \int d^{d-1} \Omega_K \simeq  K_F^{d-1} \int_{-\Lambda}^\Lambda dk \int d^{d-1} \Omega_K$, where we introduce $ k = K_F - K$ and keep only the $K_F^{d-1}$ terms after Taylor expansion. Now the kinetic part of the fermions can be re-expressed as $$\begin{aligned}
\nonumber  S_c &\simeq& K_F^{d-1} \int dk_\alpha d^{d-1} \Omega_K d\epsilon \psi^{\alpha \dagger}_c \left( i \epsilon_\alpha - v_F^\alpha k_\alpha \right) \psi^\alpha_c \\
 &=& \int d\bar{k}_\alpha d^{d-1}\bar{\Omega}_K d\bar{\epsilon} \bar{\psi}_c^{\alpha \dagger} \left ( \bar{\epsilon}_\alpha - v_F^\alpha \bar{k}_\alpha \right) \bar{\psi}_c^\alpha,\end{aligned}$$ where we introduce the dimensionless couplings, $\epsilon = \Lambda \bar{\epsilon}$, $k = \Lambda \bar{k}$, $\Omega_K = \bar{\Omega}_K$, $K_F^{d-1} \Lambda^3 \psi^\dagger \psi = \bar{\psi}^\dagger \bar{\psi}$. For the action of the bosonic fields, Eq. , we perform similar transformation, $\int d^2\bm{q} d\omega (\bm{q}^2 + \omega^2) \phi_1^2(\bm{q},\omega) = \int d^2\bar{\bm{q}} d \bar{\omega} ( \bar{\bm{q}}^2 + \bar{\omega}^2 ) \bar{\phi}_1^2 (\bar{\bm{q}},\bar{\omega})$, where we define $\Lambda^5 \phi_1^2 = \bar{\phi}_1^2$. Plugging the new definition into the Kondo action, we find that at $d=2$ it takes the form $\sim \int dk d\epsilon d^2q d\omega \psi_c^{x \dagger} (k+q, \epsilon +\omega)   \psi_c^y (k,\omega)  \omega  \phi_1 (q, \omega)$, which leads to $$\begin{aligned}
\frac{S_K}{S_c} \propto \frac{\sqrt{\Lambda}}{K_F},\end{aligned}$$ where we can see in the limit $\Lambda/K_F \rightarrow 0$, *i.e.*, the Fermi momentum is much larger than the thin-shell momentum cut-off near the Fermi surface, the Kondo coupling is heavily suppressed. Therefore, the Kondo vertex is associated with positive powers of $\sqrt{\Lambda}/K_F$ which is vanishingly small. As the number of powers of Kondo couplings increases, so does the suppression factor, and, therefore, all higher-order terms are suppressed, which means that the scaling result at tree-level RG analysis is exact. The Kondo coupling is indeed exactly *marginal*.

This exact marginality implies that the Kondo coupling does not flow to strong coupling. In other words, in the AFQ phase, the local moments do not form a multipolar Kondo singlet with the conduction electrons. Thus, the RG fixed point in the parameter regime we consider, namely weak Kondo coupling in the presence of an AFQ, shows the physics of Kondo destruction.

*Implications for the quantum phases and their transitions in heavy fermion metals*— The Kondo destruction we have shown, when the multipolar RKKY interactions dominate over the corresponding Kondo interactions, has a clear physical picture. In the AFQ order, the local degrees of freedom are strongly coupled with each other and become manifested as three quadrupolar Goldstone modes at low energies. Because these are collective bosonic modes, they can scatter the conduction electrons, but cannot form an entangled Kondo singlet with the latter. By contrast, it is well-known that when the Kondo interactions dominate over the RKKY interactions, they are marginally relevant and flow towards strong coupling, thereby yielding a Kondo entangled state ; physically, the local degrees of freedom will be able to lower the ground state energy of the system by binding with the conduction electrons into a singlet state. Calculations on the dynamical competition between RKKY and Kondo interactions from the Kondo-dominated side in Ref.  led to the proposal for two stages of Kondo destructions. Our asymptotically exact results from the opposite end shows that multipolar Kondo destruction does take place on the ordered side. As such, our results help establish a robust theoretical foundation for the notion of sequential Kondo destruction [@Martelli2017].

Our findings set the stage for detailed studies of heavy fermion materials with both spin and orbital moments in their ground state. The simplest case arises in Ce-based systems of cubic point symmetry if the $\Gamma_8$ quartet is the ground state of the $^2\rm{F}_{5/2}$ multiplet [@Pas14.1]. Examples where a continuous phase transition to a state with AFQ order was observed are, in addition to the aforementioned Ce$_3$Pd$_{20}$Si$_6$ [@Por16.1], CeB$_6$ [@Nak01.1] and, tentatively, CeTe [@Kaw11.1] and CeCoSi [@Tan18.1] under pressure.

*Conclusion*– We have studied a spin-$1$ Kondo lattice model with co-existing spin and quadrupolar local moments and used density matrix renormalization group analysis to firmly demonstrate the presence of an antiferroquadrupolar order. We have derived a non-linear sigma model description of the antiferroquadrupolar order and, based on a renormalization-group analysis, found that the Kondo couplings are exactly marginal in this phase. Our results help provide a robust theoretical foundation for the recently advanced notion of sequential localization in multipolar Kondo lattice systems  [@Martelli2017]. Our findings point to a growing list of heavy fermion metals with multipolar degrees of freedom as a new setting for the exploration towards a universal description of beyond-Landau quantum criticality and strange metal physics. In general, they illustrate how the interplay between entwined degrees of freedom can give rise to novel phases and unusual excitations, a theme that is centrally important to a broad range of strongly correlated systems.

*Acknowledgement*– The work at Rice was in part supported by the NSF (DMR-1611392), the Robert A. Welch Foundation (C-1411), the ARO (W911NF-14-1-0525), the Big-Data Private-Cloud Research Cyberinfrastructure MRI Award funded by NSF (CNS-1338099), and an IBM Shared University Research (SUR) Award. H.-H.L. has been supported by a Smalley Postdoctoral Fellowship at the Rice Center for Quantum Materials. SP acknowledges financial support from the Austrian Science Fund (project P29296-N27). Q.S. acknowledges the support of ICAM and a QuantEmX grant from the Gordon and Betty Moore Foundation through Grant No. GBMF5305, the hospitality of University of California at Berkeley and of the Aspen Center for Physics, which is supported by NSF grant No. PHY-1607611, and the hospitality and support by a Ulam Scholarship from the Center for Nonlinear Studies at Los Alamos National Laboratory.
