# Breakdown of order-fractionalization in the CPT model

Aaditya Panigrahi  Center for Materials Theory, Department of Physics and
Astronomy, Rutgers University, 136 Frelinghuysen Rd., Piscataway, NJ
08854-8019, USA  Alexei Tsvelik  Division of Condensed Matter Physics and
Materials Science, Brookhaven National Laboratory, Upton, NY 11973-5000, USA
Piers Coleman  Center for Materials Theory, Department of Physics and
Astronomy, Rutgers University, 136 Frelinghuysen Rd., Piscataway, NJ
08854-8019, USA  Department of Physics, Royal Holloway, University of London,
Egham, Surrey TW20 0EX, UK.

###### Abstract

We present an analysis of the half-filled CPT model, an analytically tractable
Kondo lattice model with Yao-Lee spin-spin interactions on a 3D hyperoctagon
lattice, proposed by Coleman, Panigrahi, and Tsvelik. Previous studies have
established that the CPT model exhibits odd-frequency triplet
superconductivity and order fractionalization. Through asymptotic analyses in
the small $J$ and large $J$ Kondo coupling limits, we identify a quantum
critical point at $J_{c}$, marking a transition from a superconductor to a
Kondo insulator. By estimating the vison gap energy to account for thermal
gauge fluctuations, we determine the energy scales governing the thermal
breakdown of order fractionalization. Moreover, at large $J$ the Kondo
insulator undergoes orbital decoupling, leading to the formation of a
decoupled Kitaev orbital liquid. These findings and analogies with the
$\mathbb{Z}_{2}$-gauged $XY$ model lead us to propose a tentative phase
diagram for the CPT model at half-filling.

###### pacs:

PACS TODO

##  I Introduction

In heavy fermion materials, the coherent scattering of conduction electrons
off a lattice of local moments produces a wide variety of emergent behavior.
These range from heavy fermion metals and superconductors [1, 2, 3, 4, 5, 6]
to topological Kondo Insulators [7, 8]. Broader classes of flat band systems,
such as Moiré materials can also be modeled as heavy fermion systems[9, 10,
11].

Figure 1: (a) The hyperoctagon lattice, a trivalent structure consisting of a
spiral of four atoms per unit cell embedded in a body centered cubic lattice.
(b) A Yao-Lee spin liquid on this lattice results in a flux-free gauge
configuration, giving rise to a Majorana Fermi surface centered at the
$P=(\pi,\pi,\pi)$ point in the Brillouin zone. (c) Schematic showing the Kondo
coupling between the conduction electrons and the local moments. (d) Cartoon
illustrating the ground-state of the large $J$ region of the CPT model,
showing in pink, the formation of Kondo singlets between conduction electrons
and spins and a decoupled Kitaev orbital liquid (KOL) .

An insightful approach to heavy fermions is to consider them as the Higgs
phase[12, 13, 14, 15, 16] of an underlying spin-liquid. From this perspective,
alternate patterns of spin fractionalization may drive new kinds of physics,
such as pairing beyond the BCS paradigm. One way to explore this idea is to
study Kondo lattices with a pre-formed spin liquid. Attempts have focused on
understanding Kondo lattice models, where the underlying spins interact via
Kitaev interactions[17, 6, 5]. Unfortunately, Kitaev-Kondo models lose their
exact solvability as a result of the Kondo term: to avoid this difficulty,
Coleman, Panigrahi, and Tsvelik (CPT) have recently proposed a three-
dimensional Kondo lattice model, [4, 3](Fig. 1) in which a Yao-Lee orbital-
spin interaction[18] restores the solvability of the Kondo lattice at half
filling.

The Hamiltonian for the CPT model on a hyperoctagon lattice (Fig 1),

$H_{CPT}=H_{c}+H_{YL}+H_{K}$ (1)

has three components:

$\displaystyle H_{c}=$ $\displaystyle-t\sum_{\left\langle
i,j\right\rangle}(c^{\dagger}_{i\sigma}c_{j\sigma}+{\rm
H.c.})-\mu\sum_{i}c^{\dagger}_{i\sigma}c_{i\sigma},$ (2) $\displaystyle
H_{YL}=$ $\displaystyle K/2\sum_{\left\langle
i,j\right\rangle}\lambda^{\alpha_{ij}}_{i}\lambda^{\alpha_{ij}}_{j}(\vec{S}_{i}\cdot\vec{S}_{j}),$
(3) $\displaystyle H_{K}=$ $\displaystyle
J\sum_{i}(c^{\dagger}_{i}\vec{\sigma}c_{i})\cdot\vec{S}_{i}$ (4)

Here $\langle i,j\rangle$ are neighboring sites on the hyper-octagonal
lattice[19], a trivalent body centered cubic (BCC) crystal with four atoms per
primitive unit cell, coiled around a helix to form alternating square and
octagonal spirals(Fig. 1b). Each site supports both electrons $c_{i\sigma}$,
localized spins $\vec{S}_{i}$ and localized orbital $\vec{\lambda}_{i}$
degrees of freedom. The conduction term $H_{c}$ describes electrons hopping
between nearest neighboring sites. The Kondo interaction $H_{K}$
antiferromagnetically couples the conduction electrons to spins $\vec{S}_{i}$
at each site. Finally, the Yao-Lee term couples the orbitals
$\vec{\lambda}_{i}$ via a Kitaev-like anisotropic interaction, “decorated” by
a Heisenberg coupling between the nearest neighbor spins $\vec{S}_{i}$. The
$\alpha_{ij}$ are the Ising coupling of orbital components along the
$\alpha_{ij}=x,y,z$ bond directions(Fig. 1c).

The anisotropic Ising coupling between orbitals induces Majorana
fractionalization [18] of spins
$\vec{S}_{j}=-(\frac{i}{2})\vec{\chi}_{j}\times\vec{\chi}_{j}$ and orbitals
$\vec{\lambda}_{i}=i\vec{b}_{j}\times\vec{b}_{j}$. In the physical Hilbert
space, where
$\sigma^{a}_{j}\lambda^{\alpha}_{j}=2i\chi^{a}_{j}b^{\alpha}_{j}$, the
fractionalized form of the Yao-Lee Hamiltonian $H_{YL}$ (3) is

$H_{YL}=K\sum_{\left\langle
i,j\right\rangle}\hat{u}_{ij}(i\vec{\chi}_{i}\cdot\vec{\chi}_{j}).$ (5)

Here, $\hat{u}_{ij}=ib^{\alpha_{ij}}_{i}b^{\alpha_{ij}}_{j}$ are the static
$\mathbb{Z}_{2}$ gauge fields, (i.e. $[H_{YL},\hat{u}_{ij}]=0$).

In three dimensions, $\mathbb{Z}_{2}$ gauge theories undergo a finite
temperature Ising phase transition at $T_{c1}$, into a deconfined phase, in
which the visons (plaquettes with a $\pi$ flux) are linearly confined. In the
Yao Lee model on a hyper-octagonal lattice, $T_{c1}\sim 0.036K$[20, 21, 4],
leading to a fractionalization of the spins into majoranas at lower
temperatures $T<T_{c1}$.

Figure 2: Proposed phase diagram of the half-filled CPT model half-filling. A
finite temperature Fermi liquid (FL) with a small Fermi surface(FL) develops
at small $J$, with a cross-over into a Kondo insulating phaseKI at large $J$.
For $T<T_{c1}$ and small $J$, the spins fractionalize into Majorana
($\vec{\chi}$) fermions, forming a 3D Yao Lee spin liquid. A logarithmic
divergence in the pairing susceptibility at small positive $J$ gives rise to a
charge $1e$, $S=1/2$ electron-majorana condensate (1eSC), with a transition
temperature comparable to the single-ion Kondo temperature $T_{K}$. Once
$T_{K}$ exceeds $T_{c1}$, $T_{c1}$ drops to zero at $J=J_{c}$, giving rise to
a superconductor-insulator quantum critical point (QCP), forming a Kondo
insulator (KI) where the spins fractionalize as Dirac ($f$) fermions. For
$T<T^{\prime}_{c1}$, the Kondo insulator co-exists with a decoupled Kitaev
orbital-liquid (KOL).

In this paper, we examine the phase diagram (Fig. 2) of the CPT model at half-
filling. Central to the phase diagram are the two characteristic energy
scales: the Kondo temperature $T_{K}$ and the Ising transition temperature
$T_{c1}$. Below the Kondo temperature $(T<T_{K})$ the conduction electrons
screen the spins, while below the Ising temperature $(T<T_{c1})$ the spins
fractionalize $\vec{S}_{j}=-\frac{i}{2}\vec{\chi}_{j}\times\vec{\chi}_{j}$
into $\vec{\chi}_{j}$ majoranas. The interplay between spin fractionalization
and Kondo screening determines the phase diagram.

Our proposed phase diagram for the half-filled CPT model (CPTM), shown in Fig.
2, is based on the analytical tractability of the extreme limits of small and
large $J$ Kondo coupling, where the model model exhibits superconductivity and
Kondo insulating behavior, respectively. One or more quantum phase transitions
at an intermediate coupling $J_{c}$ must therefore separate these distinct
phases.

Applying Ockham’s razor, we propose a single quantum critical point (QCP)
defining the superconductor-insulator transition. While we cannot entirely
rule out an alternative first-order quantum phase transition between the two
phases, this scenario would necessitate a finite temperature critical
endpoint, with a thermal rather than a quantum phase transition governing the
changing pattern of spin fractionalization.

The Kondo temperature $T_{K}(J)$ and the Ising transition temperatures
$T_{c1}(J)$ and $T_{c1}^{\prime}(J)$ divide the phase diagram (Fig. 2) into
five phases, which we summarize below:

  1. 1.

Fermi Liquid (FL) ($5t\gg T>T_{c1}$ and $T>T_{K}$): At temperatures much
higher than $T_{K}$ and $T_{c1}$, but much lower than the conduction band-
width $W\sim 5t$ the CPT model forms a conduction Fermi liquid, weakly coupled
to its embedded spins and orbitals. The spins and orbitals will exhibit
paramagnetic behavior characterized by Curie-Weiss susceptibility.

  2. 2.

Fermi Liquid* (FL∗) ($T<T_{c1}$ and $T>T_{K}$): Below the Ising transition
temperature $(T<T_{c1})$ the spins fractionalize into Majorana fermions
$\vec{S}_{j}=-\frac{i}{2}\vec{\chi}_{j}\times\vec{\chi}_{j}$, forming a Yao-
Lee spin liquid with a Majorana Fermi surface. Since the system is above the
Kondo temperature ($T>T_{K}$), the Yao-Lee spin liquid remains unscreened by
the conduction sea. Thus, in this phase, the CPT model consists of a decoupled
Fermi liquid and a Yao Lee spin liquid, constituting a Fermi liquid* phase.

  3. 3.

1e Superconductor ($1e$ SC) ($T<T_{c1}$ and $T<T_{K}$): In this phase, spins
undergo Majorana fractionalization and experience Kondo screening by
conduction electrons. The precise nesting between the Majorana Fermi surface
of the Yao-Lee spin liquid and the conduction sea at half-filling results in a
Peierls-like logarithmic instability between the electron and Majorana Fermi
surface for an infinitesimal Kondo coupling [4].

This results in an odd-frequency triplet superconductor characterized by a
fractionalized charge $e$-spinor order in the small $J$-limit [22, 3, 4]. The
system exhibits a neutral Fermi surface due to an imbalance between $4$
conduction Majoranas and $3$ spin liquid Majoranas. For more details see
section II.

  4. 4.

Kondo Insulator (KI) ($T_{K}>T>T_{c1},T^{\prime}_{c1}$ Below the Kondo
temperature $T_{K}$, but at temperatures larger than the Ising temperatures
$T_{c1},T^{\prime}_{c1}$, the spins are screened by the conduction sea,
forming a Kondo insulator, while the decoupled orbitals are unfractionalized.

  5. 5.

Kitaev Orbital Liquid (KOL) The ground state properties of the Kondo insulator
at large $J$ can be determined by carrying out a strong coupling expansion in
$1/J$. In this limit, the ground-state is a product ground state of local
Kondo singlets ($(\uparrow\Downarrow-\downarrow\Uparrow)_{j}$), co-existing
with a degenerate manifold of decoupled orbitals $\lambda^{\alpha}_{j}$

$|\Psi\rangle=\prod_{j}(\uparrow\Downarrow-\downarrow\Uparrow)_{j}|\\{\lambda_{j}\\}\rangle.$
(6)

The large spin-gap of order $J$ in this phase allows a treatment of the
Heisenberg bond operator
$\vec{S}_{i}\cdot\vec{S}_{j}\sim\langle\vec{S}_{i}\cdot\vec{S}_{j}\rangle\sim
O(t/J)$. as a static variable. Through the Yao-Lee coupling, this lifts the
orbital degeneracy, giving rise to a decoupled 3D Kitaev orbital liquid
(orbital analog of Kitaev spin liquid) below the second Ising transition
temperature $T^{\prime}_{c1}$. For more details see section III .

This paper is structured as follows: in section II we summarize the small-$J$
limit of the CPT model; in section III we present the large-$J$ limit, while
in the final section IV we discuss the intermediate coupling behavior of the
CPT model.

##  II Small-$J$ limit of the CPT model

For completeness, here we summarize the results of our previous work[4] on the
half-filled CPT model in the analytically tractable small $J$ limit.

We recall that the CPT model is a Kondo lattice model with a quartet state at
each site, comprised of a $S=1/2$ spin and orbital degree of freedom at each
site. The Yao-Lee interaction between sites acts to establish an emergent
static $\mathbb{Z}_{2}$ field[3, 4]. At half filling, the CPT model develops
nested electron and Majorana Fermi surfaces and for infinitesimal Kondo
coupling, the system undergoes a second-order phase transition into a spinor-
ordered electron-majorana condensate[4].

Yao Lee spin-liquid: The three dimensional Yao Lee (3DYL) model on a
hyperoctagon lattice shares many of the properties of a 2D Kitaev spin liquid,
most notably, the presence of gapped $Z_{2}$ flux excitations, described by
Wilson loops - products of the gauge fields $W=\prod u_{(i,j)}=\pm 1$ around
closed ten-fold loops of the hyperoctagon lattice(where $(i,j)$ orders the
sites $i$ and $j$ along $xx,\ yy$ and $zz$ bonds so that the site furthest in
the $y$, $z$ and $x$ directions respectively, is placed first [19]). In the
spin liquid ground-state, all loops are trivial $W=1$[19]; flipping the sign
of a Wilson loop creates a flux excitation (vison), with an energy determined
as a fraction of $K$. Unlike 2D Kitaev spin liquids, the 3DYL undergoes an
Ising phase transition at $T_{c1}\sim 0.036K$ into a Higgs phase where the
elementary $\mathbb{Z}_{2}$ gauge excitations (visons) are linearly confined
[23, 24, 25, 4]. and the Majorana fields describe coherent, fractionalized
spin excitations.

Below the Ising transition temperature $T_{c1}$, the confinement of visons in
the 3D Yao Lee model allows a gauge choice $u_{i,j}=1$ [20] leading to a
translationally invariant Hamiltonian. Transforming to a momentum basis,

$\vec{\chi}_{{\bf
k},\alpha}=\frac{1}{\sqrt{N}}\sum_{j}\vec{\chi}_{j\alpha}e^{-i{\bf k}\cdot{\bf
R}_{j}},$ (7)

where ${\bf R}_{j}$ is the position of the unit-cell in the BCC lattice, $N$
is the number of primitive unit cells in the lattice and $\alpha\in[1,4]$ is
the site index within each unit cell. The ground state Hamiltonian for the
Yao-Lee spin liquid is then

$H_{YL}=K\sum_{{\bf
k}\in\mbox{\mancube}}{\vec{\chi}}^{\dagger}_{\bf{k},\alpha}h({\bf
k})_{\alpha,\beta}{\vec{\chi}}_{\bf{k}\beta},$ (8)

where $\alpha,\beta\in[1,4]$ are the site indices and

$\small
h(\bf{k})=\begin{pmatrix}0&&i&&ie^{-i\bf{k}\cdot\bf{a}_{2}}&&ie^{-i\bf{k}\cdot\bf{a}_{1}}\\\
-i&&0&&-i&&ie^{-i\bf{k}\cdot\bf{a}_{3}}\\\
-ie^{i\bf{k}\cdot\bf{a}_{3}}&&i&&0&&-i\\\
-ie^{i\bf{k}\cdot\bf{a}_{1}}&&-ie^{i\bf{k}\cdot\bf{a}_{2}}&&i&&0\end{pmatrix}$
(9)

where
${\bf{a}}_{1}=(1,0,0);{\bf{a}}_{2}=\frac{1}{2}(1,1,-1);{\bf{a}}_{3}=\frac{1}{2}(1,1,1),$
are the primitive BCC lattice vectors. Since
$\chi_{-{\bf{k}}}=\chi^{\dagger}_{{\bf{k}}}$, the momentum sum is restricted
to half the Brillouin zone, corresponding to a cube () of side length $2\pi$
centered at the $P$ point at $(\pi,\pi,\pi)$. The spectrum $E_{{\bf{k}}}\equiv
K\epsilon({\bf{k}})$, determined by ${\rm
det}[\epsilon\underline{1}-h({{\bf{k}}})]=0$, or

$\epsilon^{4}-6\epsilon^{2}-8\epsilon(s_{x}s_{y}s_{z})+[9-4(s_{x}^{2}+s_{y}^{2}+s_{z}^{2})]=0,$
(10)

(where $s_{l}\equiv\sin(k_{l}/2)$, $l=x,y,z$), contains a single  Fermi
surface where $\epsilon=0$ and $s_{x}^{2}+s_{y}^{2}+s_{z}^{2}=9/4$, centered
at $P$ [26]. (Fig. 1 b.)

Conduction electrons: The nesting between Majorana and conduction Fermi
surface becomes evident upon performing the gauge transformation
$(c_{1},c_{2},c_{3},c_{4})_{\vec{R}}\rightarrow e^{i{\bf(\pi,\pi,\pi)\cdot
R}}(c_{1},ic_{2},c_{3},-ic_{4})_{\vec{R}}$ on conduction electrons in the
unit-cell at $\vec{R}$. The resulting conduction Hamiltonian $H_{c}$ takes the
form,

$H_{c}=\sum_{{\bf k}\in{\bf BZ}}c^{\dagger}_{{\bf k},\sigma\alpha}[-t\ h({\bf
k})-\mu\mathbb{I}]_{\alpha\beta}c_{{\bf k},\sigma\beta},$ (11)

where $\alpha,\beta$ denote the site indices of the unit cell, thus sharing
the same hopping matrix $h({\bf{k}})$ (8) as the spinons in the 3D Yao Lee
spin liquid. At $\mu=0$, the conduction sea develops an electron and a hole
Fermi surface at $P$ and $-P$ respectively, which can be rewritten as four
Majorana Fermi surfaces centered at $P$. These Fermi surfaces perfectly nest
with the three Majorana Fermi surface of the spin liquid, which facilitates a
BCS-like mean-field treatment of the Kondo interaction.

Kondo interaction: To this end, we express the Kondo interaction in Majorana
spin representation as,

$H_{K}=J\sum_{j}(c^{\dagger}_{j}\vec{\sigma}c_{j})\cdot(-\frac{i}{2}\vec{\chi}_{j}\times\vec{\chi}_{j})\equiv-\frac{J}{2}\sum_{l}c^{\dagger}_{j}(\vec{\chi}\cdot\vec{\sigma})^{2}c_{j}$
(12)

Small-$J$ limit: Below the Ising transition $T_{c1}$ a Hubbard-Stratonovich
transformation of the Kondo interaction in terms of the charge-$e$ spinor
order
$V_{j}=-J\langle\vec{\chi}_{a}\cdot\vec{\sigma}c_{j}\rangle=(V_{j\uparrow},V_{j\downarrow})^{T}$,

$H_{K}=\sum_{j}\bigl{[}c^{\dagger}_{j}(\vec{\sigma}\cdot\vec{\chi}_{j})V_{j}+{\rm
H.c}\bigr{]}+\frac{2V^{\dagger}_{j}V_{j}}{J}.$ (13)

The uniform saddle point of the Hamiltonian then provides us with the mean-
field solution[3, 4] in the small $J$ limit. The mismatch in the number of
conduction (4) and spin liquid majoranas (3) in equation (13) results in one
conduction majorana remaining gapless as the fractionalized order condenses.

In the mean-field solution with uniform spinor configuration
$V_{j\sigma}=(\text{V}/\sqrt{2})z_{\sigma}$, the electronic self-energy in
momentum space is expressed as:

$\Sigma_{\mathbf{k},\omega}=(1-\mathcal{Z}\otimes\mathcal{Z}^{{}^{\dagger}})V^{2}D(\mathbf{k},\omega)$
(14)

where $\mathcal{Z}=[z,i\sigma_{y}z^{*}]^{T}$ is the spinor order in Balian-
Werthamer notation, and $D(\mathbf{k},\omega)=[\omega-Kh(\mathbf{k})]^{-1}$ is
the spin-liquid propagator. The matrix
$\mathcal{Z}\otimes\mathcal{Z}^{{}^{\dagger}}$ projects out one Majorana
component of the conduction sea, leaving three components of the conduction
sea to hybridize with the spin-liquid, gapping them out in a fashion similar
to a Kondo insulator. The component of the conduction sea projected onto the
spinor ${\cal Z}$ forms a gapless “neutral” Fermi surface.

The superconducting nature of the spinor order phase is evident when
expressing the self-energy in terms of three orthogonal d-vectors, formed from
bilinears of the spinors, [4]:

$\hat{\bf{d}}^{1}+i\hat{\bf{d}}^{2}=z^{T}(-i\sigma_{2})\mathbf{\sigma}z,\hat{\bf{d}}^{3}=z^{{}^{\dagger}}{\bf{\sigma}}z.$
(15)

Written in terms of the d-vectors, the self-energy separates into magnetic and
pairing components:

$\Sigma=\Sigma_{N}({\bf k},\omega)+\Delta({\bf
k},\omega)\tau_{+}+\Delta^{{}^{\dagger}}({\bf k},\omega)\tau_{-},$ (16)

which describe the coexistence of odd-frequency magnetism and triplet
superconducting order, where

$\Sigma_{N}({\bf k},\omega)=\frac{1}{4}(3-(\hat{\bf
d}^{3}\cdot{\bf\sigma})\tau_{3})\Sigma_{0}({\bf k},\omega)$

describes odd-frequency magnetism and

$\Delta({\bf k},\omega)=-\frac{1}{4}\bigl{[}(\hat{\bf d}^{1}+i\hat{\bf
d}^{2})\cdot{\bf\sigma}\bigr{]}\Sigma_{0}({\bf k},\omega),$

describes the triplet superconductivity. Here

$\Sigma_{0}({\bf k},\omega)=V^{2}D({\bf k},\omega)$ (17)

describes the hybridization with the spin liquid. On the Fermi surface of the
spin liquid, $\Sigma_{0}({\bf k}_{F},\omega)\sim\frac{1}{\omega}$ is an odd
function of frequency. The complex d-vector
$\hat{\bf{d}}^{1}+i\hat{\bf{d}}^{2}$ breaks the time-reversal symmetry,
representing the two component superconducting order in the small-$J$ limit.

The exactness of the small-$J$ limit is a consequence of a logarithmic
divergence in the pairing susceptibility, which results in a second-order
phase transition at a transition temperature (Fig. 2),

$T_{K}=W\exp\left(-\frac{1+K/t}{\rho J}\right).$ (18)

At temperatures $T<T_{K}$, the system forms an electron-majorana pair
condensate with a spinor order parameter. Here, small $J$, combined with the
logarithmic divergence of the spinor pair susceptibility, acts as the small
parameter for the mean-field theory. Further, this spinor order generates odd-
frequency triplet pairing within the conduction sea.

##  III Large $J$ limit of the CPT Model

The large $J$ limit of the half-filled CPT model can be solved in a strong
coupling expansion in $1/J$. In this limit the half-filled ground state is a
product state of local singlets of electrons and spins at each site (Fig. 1d),
forming a Kondo insulator i.e.,

$|\Psi\rangle=\prod_{j}(\uparrow\Downarrow-\downarrow\Uparrow)_{j}|\\{\lambda_{j}\\}\rangle.$
(19)

where the ket $|\\{\lambda_{j}\\}\rangle$ describes an arbitrary configuration
of the orbital degrees of freedom. There is a gap $\Delta_{g}=J$ between the
singlet ground-state and the triplet excited states which allows us to carry
out a strong coupling expansion in powers of $K/J$ and $t/J$.

Figure 3: Depicts the Kondo insulator phase of the CPT model. The ground state
is a product state of local electron-local spin singlets, resulting in
decoupled orbitals. Through the virtual hopping of conduction electrons
between nearest neighbors, the Kondo singlet breaks into nearest neighbor
holon-doublon virtual excitation, and the spins form a singlet. The system
reverts to its original state via back-hopping to form a product state of
Kondo singlets. This virtual excitation leads to the orbitals fractionalizing
into Majorana, resulting in a decoupled Kitaev orbital liquid.

To first order in the expansion, $\langle\vec{S}_{i}\cdot\vec{S}_{j}\rangle=0$
and the expectation value of the Yao Lee term (3) is zero. However, electron
hopping (Fig. 3) gives rise to first-order corrections to the ground state:

$\left|\Psi\right\rangle=\left|\Psi_{0}\right\rangle+\sum_{j}\frac{\left|j\right\rangle\left\langle
j\right|\hat{H}_{c}\left|\Psi_{0}\right\rangle}{E_{0}-E_{j}},$ (20)

where,

$\hat{H}_{c}=-it\sum_{\left\langle
i,j\right\rangle}(c^{\dagger}_{i\sigma}c_{j\sigma}-H.c).$ (21)

Written out explicitly,

$\small\left|\Psi\\{\lambda_{j}\\}\right\rangle=\left(1-\frac{2}{3J}\left[-it\sum_{\left\langle
i,j\right\rangle}(c^{\dagger}_{i\sigma}c_{j\sigma}-H.c)\right]\right)\left|\Psi_{0}\\{\lambda_{j}\\}\right\rangle.$
(22)

The hopping moves one electron to a neighboring site forming a doubly occupied
singlet state, leading to two unscreened local moments $\vec{S}_{i}$ and
$\vec{S}_{j}$. Since the hopping preserves the spin-singlet character of the
wavefunction, the two unscreened neighboring spins must form a singlet with
$\vec{S}_{i}\cdot\vec{S}_{j}=-\frac{3}{4}$. Thus to leading order

$\left\langle\Psi\\{\lambda_{j}\\}\right|\vec{S}_{i}\cdot\vec{S}_{j}\left|\Psi\\{\lambda_{j}\\}\right\rangle=2\times\left(-\frac{3}{4}\right)\times\left(\frac{2t}{3J}\right)^{2}=-\frac{2}{3}\frac{t^{2}}{J^{2}}.$
(23)

In the low energy degenerate manifold of orbital states, the matrix elements
of the Yao-Lee interaction are

$\displaystyle\left\langle\Psi\\{\lambda_{j}\\}\right|H_{YL}\left|\Psi\\{\lambda_{j}\\}\right\rangle$
$\displaystyle=\frac{K}{2}\sum_{\left\langle
i,j\right\rangle}\left\langle\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\right\rangle\lambda^{\alpha_{ij}}_{i}\lambda^{\alpha_{ij}}_{j}$
$\displaystyle=\frac{K^{*}}{2}\sum_{\left\langle
i,j\right\rangle}\lambda^{\alpha_{ij}}_{i}\lambda^{\alpha_{ij}}_{j},$ (24)

with a renormalized coupling

$K^{*}=-\frac{8}{3}\frac{t^{2}}{J^{2}}K.$ (25)

It follows that in the large Kondo coupling $J$ limit, the ground state of the
CPT model at half-filling is a product state of spin-singlets with a decoupled
Kitaev orbital liquid.

##  IV Discussion

We now discuss the nature of the phase diagram at intermediate coupling. At
small $J$, the localized spins fractionalizes into vector ($S=1$) majoranas,
$\vec{S}_{j}=-\frac{i}{2}\vec{\chi}_{j}\times\vec{\chi}_{j}$, while at large
$J$, they fractionalize as Dirac fermions
$\vec{S}=f^{{}^{\dagger}}_{j}\vec{\sigma}f_{j}$. As the ground-state evolves
from a small $J$ superconductor to a large $J$ Kondo insulator, we anticipate
a superconductor-insulator transition. For reasons discussed below, we expect
this transition to be continuous providing a new example of a deconfined
quantum critical point, involving a transition from a $\mathbb{Z}_{2}$ Higgs
$1e$ superconductor to a $\mathbb{Z}_{2}$ deconfined Kitaev orbital liquid.

There are a number of motivating arguments for a quantum-critical
superconducting-insulator transition. Firstly, the transition involves a
Fermi-surface transformation, from a conduction neutral Majorana Fermi surface
in the superconductor to a Majorana Fermi surface of orbital excitations: this
is reminiscent of the continuous transitions from a small to a large Fermi
surface thought to occur in heavy fermion criticality[27, 28, 29]. Secondly,
energetic arguments tell us that a transition from superconductor to orbital
Kondo insulator will occur when magnetic energy of the spin liquid
($E_{mag}\sim-K$) drops below the Kondo energy gain
($E_{Kondo}\sim\frac{V^{2}}{W}$), resulting from the partial gapping of the
conduction Fermi surface, where $W$ is the band-width of the electrons. The
critical value of $K_{c}\sim\frac{V^{2}}{W}$ where this takes place implies
$T_{K}\gg K$, leading to a small vison gap energy
($\Delta_{v}\sim\frac{K^{3}}{V^{2}}$, ( see Appendix A ). This suggests a
continuous zero-temperature phase transition via the collapse of the vison gap
at the quantum critical point.

Figure 4: (a) The vison gap $\Delta_{v}$ associated with a single bond flip
$\hat{u}_{ij}=1\rightarrow-1$. $\Delta_{v}$, is the characteristic energy
scale of $\mathbb{Z}_{2}$ gauge fluctuations in the CPT model. The figure
shows the vison gap as a function of the Kondo hybridization, calculated in
the mean-field theory of the CPT model: this scale forms the boundary of the
Ising phase transition $T_{c1}$ in the regime $T_{K}>T_{c1}$. For large $J$,
the hybridization $V\sim J$, so the large $J$ behavior of vison gap energy
allows us to estimate the Ising temperature as a function of Kondo coupling
$\Delta_{v}\sim T_{c1}\sim\frac{1}{J^{2}}$.

While the vison gap energy, $\Delta_{v}\sim T_{c1}\sim\frac{1}{J^{2}}$[25]
(Fig. 4) representing the mass gap of the $\mathbb{Z}_{2}$ gauge theory,
remains finite at all $J$ in mean-field calculations ( Appendix A), there is
good reason to believe that quantum fluctuations will suppress this gap to
zero at a finite $J$, producing a continuous superconducting-insulator
transition. We expect that the zero point energy of spin-wave fluctuations,
$E_{Qfl}\sim J^{\alpha}$ is a monotonically increasing function of $J$ (i.e.
$\alpha>0$). The spin-liquid will remain energetically stable against gauge
fluctuations provided $E_{Qfl}<\Delta_{v}$. As $J$ increases, $E_{Qfl}\sim
J^{\alpha}$ rises, while $\Delta_{v}\sim 1/J^{2}$ falls continuously, so we
expect that this inequality fails at a finite critical $J_{c}$. Beyond this
point, mean-field theory breaks down due to gauge fluctuations, and the system
will transition into a Kondo insulator with a decoupled Kitaev orbital liquid.

In addition to the spin-wave fluctuations, $\mathbb{Z}_{2}$ gauge fluctuations
will also tend to suppress the stability of the $1e$ superconductor. From our
calculations of the energy spectrum associated with a single $\mathbb{Z}_{2}$
bond flip ( Appendix A: Fig. 7), we observe the associated generation of spin
$S=1$ resonant Majorana modes within the spin liquid. When the superconducting
gap $V$ exceeds these resonant energies, they sharpen into long-lived bound-
states. The concentration of these unscreened moments associated with bond-
flips will follow an activated temperature dependence,
$N\sim\exp(\frac{-(E_{bound}+\Delta_{v})}{T})$, where $E_{bound}$ is the
Majorana bound-state energy. The appearance of unscreened triplet states
implies a weakening of the Kondo screening associated with $\mathbb{Z}_{2}$
bond-flips in the superconducting phase, ultimately making the Kondo
insulator, with its robust screening, more energetically favorable.

For these reasons, we believe that the Ising phase boundaries of the $1e$
superconductor and the Kitaev orbital liquid will merge at a quantum
multicritical point. However, to definitively establish the nature of this
quantum critical point requires further numerical investigation of the model
using Monte Carlo and DMRG (on a lower dimensional analog), etc.

We now discuss the nature of the $1e$ superconducting phase and its finite
temperature superconducting phase boundary. The $1e$ superconductor is a
$\mathbb{Z}_{2}$ Higgs phase, and to understand its broken symmetry and its
finite temperature phase diagram, we must go beyond mean-field theory. The
mean-field spinor order parameter

$V(x_{j})=\langle(\vec{\sigma}\cdot\vec{\chi}_{j})c_{j}\rangle$ (26)

carries a $\mathbb{Z}_{2}$ gauge charge, and Elitzur’s theorem guarantees that
this quantity will average to zero under the thermal $\mathbb{Z}_{2}$ gauge
fluctuations. However, this does not rule out the development of charge $2e$
composite triplet order associated with

$\vec{\Psi}_{2e}=\langle V^{T}(x)i\sigma_{2}\vec{\sigma}V(x)\rangle$ (27)

which is $\mathbb{Z}_{2}$ gauge invariant, nor does it rule out the the
possibility of gauged off-diagonal long-range order, in which the charge $1e$
spinor order parameters at sites $j$ and $i$ are linked by a product of gauge
fields between the two sites,

$\langle\hat{V}^{\dagger}(x)\hat{{\cal
P}}(x,y)\hat{V}(y)\rangle\xrightarrow[]{|x-y|\rightarrow\infty}V^{\dagger}(x)V(y)$
(28)

where $\hat{{\cal P}}(x,y)=\prod_{l}u_{(l+1,l)}$ along a path from $y$ to $x$.

Figure 5: (a) Depicts a $\mathbb{Z}_{2}$ gauge XY model where the $XY$ spin
interactions are gauged by $\mathbb{Z}_{2}$ bond variables $\hat{u}_{ij}$.
Additionally, each plaquette $P_{q}$ is formed by the product of
$\hat{u}_{ij}$ around the square, costs an energy $K$ when going from
$P_{q}=1$ to $P_{q}=-1$. (b) Shows the phase diagram of $\mathbb{Z}_{2}$
gauged XY model in 3D cubic lattice. The model has 3 phases, 1. Higgs phase:
In the Higgs phase, the spin $\vec{S}$ gains long range order and vacuum
expectation value. The dashed line represents the speculated phase boundary
between the long-range ordered $e^{i\theta}$ and $e^{2i\theta}$ Higgs phases.
2. Confined phase: The confined phase is marked by the presence of deconfined
$P_{q}=-1$, $\pi$-flux plaquettes spread through the ground state. 3.
Deconfined phase: In this phase, $\pi$ flux plaquettes are absent and
$\vec{S}_{i}$s aren’t long range ordered.

These two possibilities can be explored further by integrating out out the
fermions and then expanding the effective action to leading order in the
hybridization and $\mathbb{Z}_{2}$ bond variables. The corresponding
statistical mechanical model for the finite temperature behavior is a
$\mathbb{Z}_{2}$ gauged spinor model,

$H_{SU(2)}=-J\sum_{\langle
i,j\rangle}\biggl{[}{V}^{\dagger}_{i}u_{ij}V_{j}+{\rm
H.c}\biggr{]}-K\sum_{q}\prod_{\partial P_{q}}\hat{u}_{ij}.\small$ (29)

In the presence of magnetic anisotropy or a Zeeman splitting, we can replace
$V_{j}$ by an x-y order parameter, $V_{j}\rightarrow e^{i\theta_{j}}$ , so
that this model reduces to the 3D $\mathbb{Z}_{2}$-gauged $XY$ model (or its
2+1 dimensional quantum equivalent) [30, 31, 32], where a lattice of $XY$
rotors is gauged by $\mathbb{Z}_{2}$ field (Fig 5). The Hamiltonian for the
$\mathbb{Z}_{2}$ gauged $XY$ model is given by:

$H_{qXY}=-J\sum_{\langle
ij\rangle}{u}_{ij}\cos(\theta_{i}-\theta_{j})-K\sum_{q}\prod_{\partial
P_{q}}\hat{u}_{ij}\small$ (30)

Here, each site has an $XY$ order parameter parameterized by $\theta_{i}$, and
the bonds have a $\mathbb{Z}_{2}$ gauge field $\hat{u}_{ij}$ on them (Fig.
5a), which multiplies the matter field $\theta_{i}$s. Each flux plaquette
$P_{q}$ (“vison”) at $q$ (Fig. 5a), costs an energy $2K$ in the 3D cubic
lattice. At small $J$, this model undergoes a pure gauge transition at a
finite critical $K_{c}$ in which the visons become linearly confined.

A more nuanced analysis of the 3D gauged $XY$ model (30) is necessitated to
understand the Higgs phase of the model. The J-K phase diagram of this model
is understood in the various limits $J=0$, $K=0$, $K=\infty$ and $J=\infty$.
At $J=0$, there is a deconfinement transition at a finite $K=K_{c}$ into a
phase where visons become linearly confined. At $K=0$ there is an x-y phase
transition into a phase where $\Psi_{2e}=e^{2i\theta}$ develops long range
order, whilst at $K=\infty$, where the $\mathbb{Z}_{2}$ degrees of freedom are
quenched, there is an x-y transition into a phase where $e^{i\theta}$ develops
long range order. Finally at $J=\infty$, the ground-state manifold where
$\cos(\theta_{i}-\theta_{j})=u_{ij}$, is identical to that of the
$\mathbb{Z}_{2}$ gauged Ising model considered by Fradkin and Shenkar[33],
which has no phase boundary.

What is not known about the $\mathbb{Z}_{2}$ gauged Ising model, is whether
the Ising deconfinement transition continues into the ordered phase. If the
continuation is present, then there would be a vison-confining phase
transition between a charge $2e$ and charge $1e$ order parameter. The phase
boundary of this transition would have to continue to one of the corners of
the phase diagram, that is to $J=\infty$ and either $K=0$ or $K=\infty$ (see
Fig 5b).

In drawing our tentative phase diagram of the CPT model, we have assumed that
no such phase boundary exists. However, if it does, then the Ising transition,
$T_{c1}$ would extend inside the superconducting phase, corresponding to a
transition from a charge-$e$ fractionalized ordered phase to a charge-$2e$
composite ordered phase[22]. Future numerical analysis is required to confirm
the existence of this transition.

One of the intriguing features of the proposed phase diagram (Fig. 2) are the
two routes into the order-fractionalized phase: one originating from an FL∗
phase intertwined with a decoupled spin liquid, the other emerging directly
from a heavy Fermi liquid. These two routes are reminiscent of the Bose-
Einstein and BCS condensation pathways to a conventional superconductor. This
raises the fascinating possibility that spinor superconducting order might
develop as a novel superconducting instability of a Landau Fermi liquid.
Unlike conventional triplet order, the fractionalized spinor order exhibits
Kramer’s degeneracy and a spontaneous broken time-reversal symmetry, even in
situations where there are no conventional two-dimensional triplet
representations, such as the orthorhombic triplet superconductor UTe2. At
present, the issue of whether UTe2 spontaneously breaks time-reversal symmetry
is a controversial point. Were the ground-state of this novel material to
support broken time-reversal symmetry via a single sharp phase transition,
this would constitute evidence for a fractionalized superconducting order
parameter.

###### Acknowledgements.

Acknowledgments: This work was supported by Office of Basic Energy Sciences,
Material Sciences and Engineering Division, U.S. Department of Energy (DOE)
under Contracts No. DE-SC0012704 (AMT) and DE-FG02-99ER45790 (AP and PC ).

##  Appendix A Vison gap energy and Ising phase transition

In the CPT model [4, 3], the bond variables $\hat{u}_{ij}$ commute with Kondo
interaction and remain constants of motion. Consequently, the gauge
fluctuations in the CPT model are thermal. Thermal gauge fluctuations in the
hyper-octagonal Yao Lee spin liquid subside below the Ising critical
temperature $T_{c,Ising}$. Free majoranas are the low-energy excitations in
this regime.

In the order-fractionalized phase, the condensation of electron majorana pairs
gaps out the majorana spectrum and Higgses the $\mathbb{Z}_{2}$ and Maxwell
fields. This impacts the vison gap energy $\Delta_{v}$, (the energy cost of
one bond flip) in the fractionalized ordered phase. Given that the vison gap
energy is the characteristic energy of $\mathbb{Z}_{2}$ gauge, the Ising
critical temperature varies as the vison gap energy $T_{c1}\sim\Delta_{v}$.
Indeed, Monte-Carlo simulations [25, 21] show that in 3D Kitaev spin liquids,
the Ising critical temperature $T_{c,Ising}$, and the vison gap energy
$\Delta_{v}$ are linearly correlated. Thus estimating vison gap energy
$\Delta_{v}$ provides an estimate for the Ising critical temperature.

###  A.1 Vison gap energy of CPT Model

One can analytically determine the vison gap energy for Kitaev-like spin-
liquids [34] by evaluating the change in free energy associated with flipping
a local $\mathbb{Z}_{2}$ variable $\hat{u}_{ij}$ (Fig. 6) away from the ground
state gauge configuration. A flip in the local $\mathbb{Z}_{2}$ variable
$\hat{u}_{ij}$ acts as an impurity potential, leading to a scattering phase
shift due to the Majoranas and a change in free energy, i.e., the vison gap
energy $\Delta_{v}$. This approach is validated against vison gap energy
values obtained via Monte Carlo simulations, particularly for the $V=0$ limit.
For the isotropic case $J_{x}=J_{y}=J_{z}=K$, Monte Carlo simulations yield a
vison gap energy of $\Delta_{v}=0.09(1)K$[21] for the Kitaev spin liquid on a
hyperoctagon lattice (i.e. (10,3)a system). Our calculation aligns precisely
with this result, yielding a vison gap energy of $\Delta_{v}=0.089(5)K$ per
Majorana species $\chi^{a}$ i.e. the Yao-Lee equivalent of Kitaev spin liquid.
This methodology remains robust even when the Majoranas hybridize with
electrons, allowing us to extend the approach to the CPT model.

Figure 6: (a) Shows the plaquettes in hyperoctagon lattice that change signs
when one flips a bond variable $\hat{u}_{ij}=1\rightarrow-1$ shown in the
oval, the energy. The energy cost associated with this flip is the vison gap
energy (b)Vison gap energy as a function of hybridization calculated
numerically on a 50x50x50 unit cell lattice, with a hybridization resolution
of $\frac{\Delta V}{K}=0.01$. Consequently, we find the best fit for the vison
gap energy as a function of hybridization $V$.

To compute the vison gap energy in the CPT model, we flip the $\hat{u}_{ij}$
bond between the $2^{nd}$ and $3^{rd}$ atoms in the unit cell positioned at
the origin,

$H_{CPT+3v}=H_{CPT}-2K(i\vec{\chi}_{0,2}\cdot\vec{\chi}_{0,3}).$ (31)

$H_{CPT+3v}$ is the Hamiltonian that has a local $\mathbb{Z}_{2}$ bond flip
over this ground-state gauge configuration. This bond flip is associated with
the creation of 3 visons adjacent to the bond. Treating the bond-flip term as
an impurity potential,

$\hat{V}_{flip}=-2K(i\vec{\chi}_{0,2}\cdot\vec{\chi}_{0,3})$ (32)

allows us to calculate the associated free-energy change. The free energy of
the CPT model with $\mathbb{Z}_{2}$ bond flip is expressed in terms of its
bare-Green’s function $G_{0,CPT}$,

$\beta F=-\frac{1}{2}Tr\log[-G_{0,CPT}^{-1}+\hat{V}_{flip}]$ (33)

where $G_{0,CPT}(\omega,\vec{k})=(\omega\mathbb{I}-h_{CPT}(\vec{k}))^{-1}$ is
the Green’s function for the CPT Model in the mean-field configuration.

Since the bond-flip potential $\hat{V}_{flip}$ scatters Majorana fermions in
the Yao-Lee spin liquid. The associated free-energy change is obtained in
terms of effective majorana green’s function $G_{\vec{\chi}}$, which includes
self-energy corrections from the electron-majorana condensate, as follows

$\Delta F=\frac{1}{2\beta}Tr\log[1-\hat{V}_{flip}G_{\vec{\chi}}].$ (34)

Here, the trace is over the system and Matsubara frequencies.

The effective Majorana Green’s function in the Majorana-electron condensate is
given by,

$\large
G_{\vec{\chi}}(z,\vec{k},V)=\frac{1}{z-Kh(\vec{k})-\frac{V^{2}}{z+t\,h(\vec{k})}}$
(35)

Where $V$ is the magnitude of the spinor order, and $h(\vec{k})$ is the
$4$-band Hamiltonian given in equation (8). This effective Green’s function is
used to calculate vison gap energy $\Delta_{v}$ by calculating the scattering
phase shift of the $\mathbb{Z}_{2}$ bond-flip potential.

The scattering potential $\hat{V}_{flip}$ is local, the free-energy of
$\mathbb{Z}_{2}$ bond flip potential can be re-expressed in terms of the local
majorana Green’s function $g_{\vec{\chi}}(i\omega_{n})$

$\Delta
F=\frac{1}{2\beta}\sum_{i\omega_{n}}[\ln(1-\hat{V}_{flip}g_{\vec{\chi}}(i\omega_{n}))]$
(36)

where, the local majorana Green’s function $g_{\vec{\chi}}(i\omega_{n})$ is,

$g_{\vec{\chi}}(z)=\frac{1}{N_{c}}\sum_{k\in BZ}G_{\vec{\chi}}(z,k)$ (37)

obtained by summing the Majorana Green’s function over the Brillouin zone.

Figure 7: (a) Depicts the scattering phase shift associated with a bond flip
impurity $\hat{V}_{flip}$ for a species of spin Majorana $\chi^{a}$,
illustrating two spinor order strengths: $V=0K$ (in blue) and $V=2.9K$ (in
red). The plateau for the gapped $V=2.9K$ state at $\Delta_{v}=\pi$ signifies
the formation of an in-gap bound state.(b)Shows the change in density of
states $\delta\rho=\frac{1}{2\pi}\frac{d\delta_{v}(\omega)}{d\omega}$ (in
blue) associated with the bond flip potential $\hat{V}_{flip}$. The sharp peak
inside the gap (shown using the bulk density of state $\rho$ for $V=2.9K$ in
red) signifies the presence of an in-gap bound state for spinor order strength
of $V=2.9K$.

Upon carrying out the Matsubara frequency summation, one obtains the
expression of the free-energy change in terms of the scattering phase shift,

$\Delta
F=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\left(\frac{1}{2}-f(\omega)\right)\delta_{v}(\omega)$
(38)

where, the scattering phase shift $\delta_{v}$ is given by

$\delta_{v}(\omega)=Im{\hbox{Tr}}\log[1-\hat{V}_{flip}g_{\vec{\chi}}(z)]_{z=\omega-i\delta}.$
(39)

At zero-temperature, this free-energy change corresponds to the vison gap
energy $\Delta_{v}$ given by,

$\Delta_{v}=-K\int_{-\infty}^{0}\frac{dx}{2\pi}Im\log[\det(1-\hat{V}_{flip}g_{\vec{\chi}}(z))]_{z=x-i\delta}.$
(40)

Where the Tr is over the sites within the unit-cells. Numerically, the vison
gap energy calculation was carried out on a $50\times 50\times 50$ lattice by
carrying out discrete summation over the momentum in the Brillouin zone and
frequencies.

We compute the vison gap energy for discrete values of hybridization $V$ with
a resolution of $\frac{\Delta V}{t}=0.01$. The resulting dataset was fit to
obtain a functional form,

$\Delta_{v}(V)=r\left(\frac{1}{(sV+u)^{2}}-\frac{1}{(sV+u)^{4}}\right)$ (41)

with $s=0.273$, $u=1.040$ and $r=3.680$. This function form matches the
asymptotic behavior of the vison gap energy $\Delta_{v}\sim\frac{1}{V^{2}}$
for large $V$, which is always positive and is obtained analytically using
equation (40). Additionally, it also matches the numerical result [21] for
$V=0$ as can be seen in (Fig. 6), and consequently is a reliable fit. The
initial enhancement in the vison gap energy $\Delta_{v}$ is a consequence of
$\mathbb{Z}_{2}$ gauge theories being Higgs phases of continuum gauge theories
[35, 32], thus electron-majorana fractionalized order formation enhances the
already massive $\mathbb{Z}_{2}$ gauge fields. The reduction in the vison gap
energy $\Delta_{v}$ for large hybridization $\Delta_{v}$ is a result of
renormalization effects at large Kondo coupling $J$.

Since the vison gap energy in the CPT model asymptotically decreases as
$\Delta_{v}\sim\frac{1}{V^{2}}$ for large hybridization $V$. Thus, as the
Kondo coupling $J$ increases, the characteristic vison gap energy scale
renormalizes to smaller values as a result of Kondo screening. Additionally,
given the linear correlation between vison-gap energy and Ising critical
temperature $T_{c,Ising}$, we estimate that Ising critical temperature
$T_{c,Ising}\sim\Delta_{v}\sim\frac{1}{V^{2}}$ reduces with increasing Kondo
coupling. Beyond this Ising critical temperature scale $T_{c,Ising}$, the
thermal gauge fluctuations destroy electron-majorana fractionalized order.

The electron-majorana fractionalized order in the CPT model is expected to be
suppressed by quantum fluctuations about the mean-field theory. Such quantum
fluctuations will grow as one moves away fromthe small Kondo coupling $J$
regime, where the $J$ is the small parameter that controlling the mean-field
treatment. In the large Kondo coupling $J$ limit, the ground state is a Kondo
insulator with a decoupled orbital Kitaev spin liquid. Thus, passing from the
superconducting electron-majorana condensate at small $J$ to the Kondo
insulator phase at large $J$, the system undergoes a quantum phase transition.
Away from half-filling, this quantum phase transition is associated with the
small-to-large expansion of the neutral Fermi surface, a likely signature of a
continuous quantum phase transition.

## References

  * Petrovic  _et al._ [2001] C. Petrovic, P. G. Pagliuso, M. F. Hundley, R. Movshovich, J. L. Sarrao, J. D. Thompson, Z. Fisk, and P. Monthoux, Journal of Physics: Condensed Matter 13, L337 (2001).
  * Flint  _et al._ [2008] R. Flint, M. Dzero, and P. Coleman, Nature Physics 4, 643 (2008).
  * Tsvelik and Coleman [2022] A. M. Tsvelik and P. Coleman, Physical Review B 106, 125144 (2022).
  * Coleman  _et al._ [2022] P. Coleman, A. Panigrahi, and A. Tsvelik, Physical Review Letters 129, 177601 (2022).
  * Choi  _et al._ [2018] W. Choi, P. W. Klein, A. Rosch, and Y. B. Kim, Physical Review B 98, 155123 (2018).
  * Seifert  _et al._ [2018] U. F. P. Seifert, T. Meng, and M. Vojta, Physical Review B 97, 085118 (2018).
  * Dzero  _et al._ [2010] M. Dzero, K. Sun, V. Galitski, and P. Coleman, Physical Review Letters 104, 106408 (2010).
  * Neupane  _et al._ [2013] M. Neupane, N. Alidoust, S.-Y. Xu, T. Kondo, Y. Ishida, D. J. Kim, C. Liu, I. Belopolski, Y. J. Jo, T.-R. Chang, H.-T. Jeng, T. Durakiewicz, L. Balicas, H. Lin, A. Bansil, S. Shin, Z. Fisk, and M. Z. Hasan, Nature Communications 4, 2991 (2013).
  * Song and Bernevig [2022] Z.-D. Song and B. A. Bernevig, Physical Review Letters 129, 047601 (2022).
  * Kumar  _et al._ [2022] A. Kumar, N. C. Hu, A. H. MacDonald, and A. C. Potter, Physical Review B 106, L041116 (2022).
  * Ramires and Lado [2018] A. Ramires and J. L. Lado, Physical Review Letters 121, 146801 (2018).
  * Coleman [1983] P. Coleman, Physical Review B 28, 5255 (1983).
  * Read and Newns [1983] N. Read and D. M. Newns, Journal of Physics C: Solid State Physics 16, 3273 (1983).
  * Kogut [1979] J. B. Kogut, Reviews of Modern Physics 51, 659 (1979).
  * Senthil  _et al._ [2003] T. Senthil, M. Vojta, and S. Sachdev, Phys. Rev. Lett 90, 216403 (2003).
  * Coleman  _et al._ [2005] P. Coleman, J. B. Marston, and A. J. Schofield, Phys. Rev. B 72, 245111. 15 p (2005).
  * Kitaev [2006] A. Kitaev, Annals of Physics 321, 2 (2006).
  * Yao and Lee [2011] H. Yao and D.-H. Lee, Physical Review Letters 107, 087205 (2011).
  * Hermanns and Trebst [2014] M. Hermanns and S. Trebst, Phys. Rev. B 89, 235102 (2014).
  * Hermanns  _et al._ [2018] M. Hermanns, I. Kimchi, and J. Knolle, Annual Review of Condensed Matter Physics 9, 17 (2018).
  * O’Brien  _et al._ [2016] K. O’Brien, M. Hermanns, and S. Trebst, Physical Review B 93, 085101 (2016).
  * Komijani  _et al._ [2018] Y. Komijani, A. Toth, P. Chandra, and P. Coleman,  (2018), 10.48550/ARXIV.1811.11115, publisher: arXiv Version Number: 2.
  * Seifert  _et al._ [2020] U. F. P. Seifert, X.-Y. Dong, S. Chulliparambil, M. Vojta, H.-H. Tu, and L. Janssen, Phys. Rev. Lett. 125, 257202 (2020).
  * Mishchenko  _et al._ [2017] P. A. Mishchenko, Y. Kato, and Y. Motome, Phys. Rev. B 96, 125124 (2017).
  * Eschmann  _et al._ [2020] T. Eschmann, P. A. Mishchenko, K. O’Brien, T. A. Bojesen, Y. Kato, M. Hermanns, Y. Motome, and S. Trebst, Physical Review B 102, 075125 (2020).
  * Hermanns  _et al._ [2015] M. Hermanns, S. Trebst, and A. Rosch, Phys. Rev. Lett. 115, 177205 (2015).
  * Coleman  _et al._ [2001] P. Coleman, C. Pépin, Q. Si, and R. Ramazashvili, Journal of Physics: Condensed Matter 13, R723 (2001).
  * Senthil  _et al._ [2004] T. Senthil, M. Vojta, and S. Sachdev, Physical Review B 69, 035111 (2004).
  * 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).
  * Lammert  _et al._ [1993] P. E. Lammert, D. S. Rokhsar, and J. Toner, Physical Review Letters 70, 1650 (1993).
  * Lammert  _et al._ [1995] P. E. Lammert, D. S. Rokhsar, and J. Toner, Physical Review E 52, 1778 (1995).
  * Sachdev [2019] S. Sachdev, Reports on Progress in Physics 82, 014001 (2019).
  * Fradkin and Shenker [1979] E. Fradkin and S. H. Shenker, Physical Review D 19, 3682 (1979).
  * Panigrahi  _et al._ [2023] A. Panigrahi, P. Coleman, and A. Tsvelik, Physical Review B 108, 045151 (2023).
  * Fradkin [2021] E. Fradkin,  _Quantum field theory: an integrated approach_ (Princeton University Press, Princeton, 2021).


# Price-Aware Deep Learning for Electricity Markets

Vladimir Dvorkin, and Ferdinando Fioretto

###### Abstract

While deep learning gradually penetrates operational planning, its inherent
prediction errors may significantly affect electricity prices. This letter
examines how prediction errors propagate into electricity prices, revealing
notable pricing errors and their spatial disparity in congested power systems.
To improve fairness, we propose to embed electricity market-clearing
optimization as a deep learning layer. Differentiating through this layer
allows for balancing between prediction and pricing errors, as oppose to
minimizing prediction errors alone. This layer implicitly optimizes fairness
and controls the spatial distribution of price errors across the system. We
showcase the price-aware deep learning in the nexus of wind power forecasting
and short-term electricity market clearing.

###### Index Terms:

Deep learning, differentiable programming, electricity markets, fairness,
optimal power flow, wind power forecasting

##  I Introduction

Addressing the gap between electricity market clearing and real-time
operations, and to accelerate market-clearing subroutines, recent literature
proposes various deep learning tools for forecasting weather-dependent
generation and loads, set-points of conventional generation, and active sets
of network constraints, among other applications. Deep learning models learn
complex physical processes from an abundance of operational data and instantly
map contextual features into informed inputs for market clearing [1], and even
directly estimate market-clearing outcomes from data [2]. However well-
informed and tuned these models may be, they are still prone to prediction
errors, eventually affecting market outcomes for electricity producers and
consumers.

This letter uncovers the propagation of prediction errors of deep learning
into locational marginal prices – the ultimate indicator defining electricity
market outcomes across the network – and demonstrates spatial disparity of
these errors across the range of power system benchmark systems. To reduce the
errors and improve fairness in electricity market applications, we propose to
directly incorporate market-clearing optimization as a deep learning layer to
inform predictions on pricing errors. For clarity, we use wind power
forecasting as the deep learning task, but other prediction tasks may benefit
from the same approach. Our results show that even small errors cause spatial
disparities in congested systems, with usually a small set of buses
contributing the most to the overall disparity. While transmission investments
would ultimately resolve spatial disparities, informing predictions with
market-clearing outcomes already improves fairness.

##  II Price-Awareness for Deep Learning

Predicting wind power output, we use a dataset
$\\{(\varphi_{i},w_{i})\\}_{i=1}^{m}$ of $m$ operational records, where the
$i^{\text{th}}$ record includes a power output $w_{i}$ and the associated
feature vector $\varphi_{i}$, collecting weather and turbine data, such as
wind speeds and blade pitch angles. This dataset is used to train a deep
learning model DeepWP that maps features to wind power outputs, i.e.,
$\texttt{DeepWP}(\varphi)=\widehat{w}.$ The model is restricted to neural
networks with a fully connected feedforward architecture, and we train it to
minimize the prediction error $\lVert\widehat{w}-w\rVert.$ While standard,
this training approach is unaware of the effects of prediction errors on
electricity prices.

###  II-A From Prediction to Electricity Price Errors

To study the impact of prediction errors on electricity prices, we work with a
simplified model of electricity market clearing based on a DC optimal power
flow (OPF) problem [1, 2]. Given a wind power forecast
$\widehat{w}\in\mathbb{R}^{n}$, the DC-OPF goal is to identify the
conventional generation dispatch $p\in\mathbb{R}^{n}$ which satisfies electric
loads $d\in\mathbb{R}_{+}^{n}$ at minimum dispatch costs. Generators produce
within technical limits $\underline{p},\overline{p}\in\mathbb{R}_{+}^{n}$
incurring costs, described by a quadratic function with coefficients
$c\in\mathbb{R}_{+}^{n}$ and $C\in\mathbb{S}_{+}^{n}$. The power flows are
computed using a matrix $F\in\mathbb{R}^{e\times n}$ of power transfer
distribution factors, i.e., the flows are
$f=F(p+\widehat{w}-d)\in\mathbb{R}^{e}$ limited by
$\overline{f}\in\mathbb{R}_{+}^{e}$. The DC-OPF then takes the form:

$\displaystyle\underset{{\underline{p}\leqslant
p\leqslant\overline{p}}}{\text{minimize}}\quad$ $\displaystyle
p^{\top}Cp+c^{\top}p$ (1a) subject to
$\displaystyle\mathbb{1}^{\top}(p+\widehat{w}-d)=0\;\colon\widehat{\lambda}_{b},$
(1b)
$\displaystyle|F(p+\widehat{w}-d)|\leqslant\overline{f}\;\colon\widehat{\lambda}_{\overline{f}},\widehat{\lambda}_{\underline{f}},$
(1c)

which minimizes the total generation cost, while respecting power balance and
the limits on power generation and transmission.

The locational marginal prices (LMPs), induced on a particular wind power
forecast $\widehat{w},$ are derived from the optimal dual variables accociated
with power balance and power flow constraints:

$\displaystyle\pi(\widehat{w})=\widehat{\lambda}_{b}\cdot\mathbb{1}-F^{\top}(\widehat{\lambda}_{\overline{f}}-\widehat{\lambda}_{\underline{f}})\in\mathbb{R}^{n},$
(2)

where the first component is a system-wide price and the second component
corrects this price to account for network congestion.

######  Remark 1 (Uniqueness).

If convex problem (1) is feasible for some forecast $\widehat{w}$, then its
primal and dual solutions are unique thanks to strictly monotone objective
function (1) [3]. Hence, electricity price $\pi(\widehat{w})$ is also unique
with respect to forecast $\widehat{w}$.

Consequently, the LMP error, which is defined as

$\displaystyle\delta\pi=\pi(\widehat{w})-\pi(w)\in\mathbb{R}^{n},$ (3)

is also unique with respect to wind power forecast $\widehat{w}$.

Two key observations arise from analyzing Eq. (2). The first relates to
disparities arising from network congestion.

######  Observation 1 (Spatial disparity).

In congested networks, for which
$\mathbb{1}^{\top}(\widehat{\lambda}_{\overline{f}}+\widehat{\lambda}_{\underline{f}})>0$,
the price error at bus $i$ is proportional to the $i^{\text{th}}$ column of
matrix $F$ of power transfer distribution factors.

Hence, the same prediction error has a disparate effect on electricity prices,
depending on location. This effect under DeepWP predictions is illustrated in
Fig. 2 (left plane), where LMP errors vary in the range from $-4.2$ $/MWh to
$+0.9$ $/MWh on average.

The second observation follows from the first one and identifies a unique
system bus, reference bus $r$, with the smallest price error.

######  Observation 2 (Reference bus).

Since the $r^{th}$ column of $F$ is all zeros, the price error at the
reference bus only includes the error of the system-wide component in (2).

The spatial disparities across system buses is measured using the notion of
$\alpha-$fairness [4]:

$\displaystyle\alpha=\underset{i\in
1,\dots,n}{\text{max}}\lVert\mathbb{E}[\lVert\delta\pi_{i}\rVert]-\mathbb{E}[\lVert\delta\pi_{r}\rVert]\rVert.$
(4)

where the expectation is with respect to the dataset distribution. Parameter
$\alpha$ is called the fairness bound, with smaller values denoting stronger
fairness. Although the original definition in [4] takes the maximum across all
pairs of individuals (buses), we use the fact that the reference bus $r$ has
the smallest error.

###  II-B Electricity Market Clearing as a Neural Network Layer

To minimize the impact of prediction errors on electricity prices, we develop
a new deep learning architecture DeepWP+, depicted in Fig. 1. Relative to
DeepWP, the new architecture includes a market-clearing optimization layer
which computes electricity prices in response to wind power predictions. The
prices then enter the loss function to minimize a linear combination of wind
power prediction and electricity price errors.

Directly acting on problem (1), however, is computationally challenging: it
requires solving problem (1) and then quering the dual solution – operation
which does not admit an explicit gradient. We propose to work with the dual
optimization instead. For brevity, we rewrite problem (1) in a compact form:

$\displaystyle\underset{{x}}{\text{minimize}}\quad$ $\displaystyle
p^{\top}Cp+c^{\top}p$ (5a) subject to $\displaystyle Ap\geqslant
b(\widehat{w})\;\colon\widehat{\lambda},$ (5b)

where the feasible region is parameterized by wind power forecast
$\widehat{w}$. The corresponding dual problem takes the form:

$\displaystyle\underset{{\lambda\geqslant\mathbb{0}}}{\text{maximize}}\quad$
$\displaystyle q(\widehat{w})^{\top}\lambda-\lambda^{\top}Q\lambda,$ (6)

where $q(\widehat{w})=AC^{-1}c+b(\widehat{w})$ and $Q=AC^{-1}A^{\top}$. Since
the two problems are strictly convex and concave, respectively, the strong
duality holds, such that the LMPs can be extracted from the dual problem using
expression (2). Another virtue of the dual problem (6) is that it only
includes the non-negativity constraints, and it is thus simpler to solve than
its primal counterpart (5).

$\widehat{w}$features $\varphi$ReLUReLUReLUforecastdual of DC-OPF
(6)$\widehat{w}$$\pi(\widehat{w})$$\lVert\widehat{w}-w\rVert+\lVert\pi(\widehat{w})-\pi(w)\rVert$DeepWPDeepWP+
Figure 1: The standard DeepWP and proposed DeepWP+ learning architectures.

The market clearing as an optimization layer is implemented using DiffOpt.jl –
a library for differentiating through the solution of optimization problems in
Julia [5].

##  III Numerical Results and Discussion

For numerical experiments, we use power system benchmark systems from [6]. In
each system, we install one wind farm and uniformly scale transmission
capacity to provoke congestion, as shown in Table I. To ensure fair
compassion, we use the same wind power forecasting
data111https://www.kaggle.com/datasets/theforcecoder/wind-power-forecasting
across all systems. The data includes wind power output as a function of wind
speed, wind direction and blade pitch angle features. We independently sample
1,000 scenarios for training and testing.

TABLE I: Wind Power prediction and LMP errors under conventional (DeepWP) and price-aware (DeepWP+) deep learning models case | wind power data | DeepWP | DeepWP+
---|---|---|---
bus | capacity | $\overline{f}$-scale | RMSE$(\widehat{w})$ | RMSE$(\widehat{\pi})$ | CVaR$(\widehat{\pi})$ | $\alpha-$value | RMSE$(\widehat{w})$ | RMSE$(\widehat{\pi})$ | CVaR$(\widehat{\pi})$ |  $\alpha-$value
MW | [p.u.] | MWh | $/MWh | $/MWh | $/MWh | MWh | gain | $/MWh | gain | $/MWh | gain | $/MWh | gain
14_ieee | 14 | 100 | 1.00 | 0.35 | 0.62 | 1.52 | 0 | 0.35 | $+0.6\%$ | 0.61 | $\mathbf{-0.6\%}$ | 1.50 | $-0.8\%$ | 0 | —
57_ieee | 38 | 600 | 0.60 | 2.31 | 11.03 | 34.64 | 32.08 | 2.60 | $+11.2\%$ | 10.72 | $\mathbf{-2.9\%}$ | 33.59 | $-3.1\%$ | 30.92 | $-3.8\%$
24_ieee | 15 | 1,000 | 0.75 | 4.08 | 8.62 | 37.70 | 27.48 | 4.51 | $+9.6\%$ | 8.33 | $\mathbf{-3.5\%}$ | 36.35 | $-3.7\%$ | 26.26 | $-4.6\%$
39_epri | 6 | 1,500 | 0.70 | 5.94 | 11.15 | 31.21 | 17.53 | 6.43 | $+7.6\%$ | 10.19 | $\mathbf{-9.4\%}$ | 28.02 | $-11.4\%$ | 15.84 | $-10.7\%$
73_ieee | 41 | 1,000 | 0.80 | 4.02 | 5.12 | 16.21 | 32.83 | 5.51 | $+26.9\%$ | 4.24 | $\mathbf{-20.8\%}$ | 13.41 | $-20.9\%$ | 26.63 | $-23.3\%$
118_ieee | 37 | 500 | 0.75 | 2.29 | 3.59 | 11.32 | 17.91 | 2.60 | $+12.1\%$ | 2.88 | $\mathbf{-24.7\%}$ | 9.06 | $-25.0\%$ | 14.09 | $-27.2\%$
Figure 2: Projection of conventional DeepWP and price-aware DeepWP+ wind power
predictions on locational marginal price errors in the IEEE 118-bus system.

The DeepWP architecture consists of four hidden layers with 30 neurons each,
all using ReLU as activation functions. The training consists of three stages,
all using ADAM optimizer: pre-training with 500 epochs and learning rate
$1e$–$4$, then 1000 epochs with learning rate $5e$–$5$, and final 100 epochs
with learning rate $5e$–$6$. The training of DeepWP+ starts from the
$501^{\text{th}}$ epoch using the parameters of the pre-trained DeepWP model
obtained at the first stage. As neural network parameters are initialized at
random, we report the average results across 100 trained models.

We discuss the results using the root mean square error of wind power
prediction, denoted $\text{RMSE}(\widehat{w})$, that of prices, denoted
$\text{RMSE}(\widehat{\pi})$, the root mean square error across 10% of the
worst-case scenarios, denoted $\text{CVaR}(\widehat{\pi})$, and
$\alpha-$fairness bound.

Table I provides a statistical summary for different systems. Under prediction
errors of the price-agnostic DeepWP model, LMP errors vary between 0.62 and
11.15 $/MWh This range, obtained using the same wind power dataset, is
explained by at least two factors. First, all systems have unique cost
structures: the average LMP spans the rage from 42.8 (14_ieee) to 149.3
(39_epri) $/MWh, for the median wind power output. Second, congested networks
prone to larger errors. This is particularly evident from the 14_ieee system,
where no wind power output scenario causes congestion. The zero fairness bound
$\alpha$ for this system also supports this observation.

The application of the price-aware DeepWP+ model, on the other hand,
demonstrates the reduction of LMP errors, varying from 0.6 to 24.7% relative
to the DeepWP model. Even larger reductions are observed across 10% of the
worst-case wind scenarios, measured by $\text{CVaR}(\widehat{\pi})$. For
example, in the 39_epri case, while the average LMP error is reduced by 0.96
$/MWh, the worst-case error decreases by 3.19 $/MWh. Although the focus of
DeepWP+ is on price errors, it implicitly minimized the disparity of LMP
errors, as shown by the reduced fairness bound $\alpha$ across all congested
systems. Figure 2 depicts this effect for the most geographically distributed
118_ieee system, where the fairness bound $\alpha$ is reduced by 27.2% thanks
to price-aware predictions. Observe, while the majority of buses demonstrate
near zero price errors, it is a small set of buses contributing to
$\alpha-$fairness statistic the most. Importantly, the benefits of DeepWP+
predictions do come at the expense of increasing prediction errors. However,
all DeepWP+ predictions are feasible and yield competitive and fairer – in the
$\alpha-$fairness sense – electricity prices.

Figure 3: Prediction (top) and locational marginal price (bottom) errors in the 39-bus EPRI system during the training of DeepWP and DeepWP+ models. TABLE II: Average CPU time per training epoch of DeepWP+ architecture case | 14_ieee | 24_ieee | 39_epri | 57_ieee | 73_ieee | 118_ieee
---|---|---|---|---|---|---
time [min.] | 0.06 | 0.26 | 0.22 | 0.48 | 2.03 | 6.28

The price error reduction and fairness improvements are due to the bias which
DeepWP+ introduces during the training procedure. Figure 3 visualizes how the
bias is introduced starting from the $501^{\text{th}}$ epoch, which diverts
the training towards the desired result, i.e., to minimal prediction or
minimal price error. Finally, we remark that incorporating the electricity
market as an optimization layer in deep learning increases the computational
burden of the training procedure. Table II reports CPU times per training
epoch, which tends to increase in the size of the system.

Overall, our results show that embedding market clearing as a deep learning
layer informs predictions on market outcomes and improves algorithmic fairness
in electricity markets.

## References

  * [1] R. Ferrando _et al._ , “A physics-informed machine learning for electricity markets: A NYISO case study,” _arXiv preprint arXiv:2304.00062_ , 2023.
  * [2] S. Liu _et al._ , “Topology-aware graph neural networks for learning feasible and adaptive AC-OPF solutions,” _IEEE Trans. Power Syst._ , 2022\.
  * [3] V. Dvorkin Jr _et al._ , “Electricity market equilibrium under information asymmetry,” _Oper. Res. Lett._ , vol. 47, no. 6, pp. 521–526, 2019.
  * [4] C. Tran _et al._ , “Decision making with differential privacy under the fairness lens,” in _International Joint Conference on Artificial Intelligence_ , 2021.
  * [5] A. Sharma _et al._ , “Flexible differentiable optimization via model transformations,” _arXiv preprint arXiv:2206.06135_ , 2022.
  * [6] C. Coffrin _et al._ , “Powermodels. jl: An open-source framework for exploring power flow formulations,” in _Power Systems Computation Conference_ , 2018.