# Hilbert curve vs Hilbert space:
exploiting fractal 2D covering to increase tensor network efficiency

Giovanni Cataldi  Dipartimento di Fisica e Astronomia “G. Galilei”, Università
di Padova, I-35131 Padova, Italy.  Padua Quantum Technologies Research Center,
Università degli Studi di Padova.  Istituto Nazionale di Fisica Nucleare
(INFN), Sezione di Padova, I-35131 Padova, Italy.  Ashkan Abedi  Scuola
Normale Superiore, I-56127 Pisa, Italy.  Giuseppe Magnifico  Dipartimento di
Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova, Italy.
Padua Quantum Technologies Research Center, Università degli Studi di Padova.
Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131
Padova, Italy.  Simone Notarnicola  Dipartimento di Fisica e Astronomia “G.
Galilei”, Università di Padova, I-35131 Padova, Italy.  Padua Quantum
Technologies Research Center, Università degli Studi di Padova.  Istituto
Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy.
Nicola Dalla Pozza  Scuola Normale Superiore, I-56127 Pisa, Italy.  Vittorio
Giovannetti  NEST, Scuola Normale Superiore and Istituto Nanoscienze-CNR,
I-56127 Pisa, Italy.  Simone Montangero  Dipartimento di Fisica e Astronomia
“G. Galilei”, Università di Padova, I-35131 Padova, Italy.  Padua Quantum
Technologies Research Center, Università degli Studi di Padova.  Istituto
Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy.

(August 27, 2024)

###### Abstract

We present a novel mapping for studying 2D many-body quantum systems by
solving an effective, one-dimensional long-range model in place of the
original two-dimensional short-range one. In particular, we address the
problem of choosing an efficient mapping from the 2D lattice to a 1D chain
that preserves at most the locality of interactions within the TN structure.
By using Matrix Product States (MPS) and Tree Tensor Network (TTN) algorithms,
we compute the ground state of the 2D quantum Ising model in transverse field
with lattice size up to $64\times 64$, comparing the results obtained from
different mappings based on two space-filling curves, the snake curve and the
Hilbert curve. We show that the locality-preserving properties of the Hilbert
curve leads to a clear improvement of numerical precision, especially for
large sizes, and turns out to provide the best performances for the simulation
of 2D lattice systems via 1D TN structures.

Unveiling the static and dynamical properties of strongly correlated quantum
many-body systems remains one of the fundamental goals of the current research
in condensed-matter physics. Quantum systems with a large number of
interacting microscopic constituents can generate collective states that have
no classical counterpart and are extremely interesting both from a theoretical
point of view as well as for possible experimental applications [1, 2, 3].
Paradigmatic examples include quantum Hall states [4], superconductors [5],
topological insulators [6] and quantum spin liquids [7]. Such collective
quantum phenomena play a central role not only in condensed-matter but also in
high-energy physics, quantum chemistry, atomic physics as well as in quantum
technologies where it is of paramount importance to understand and control the
quantum states of many-body systems for quantum simulation and computation [8,
9, 10].

(a)

(b)

(c)

(d)

Figure 1: Pictorial representation of the Snake (a) and the Hilbert curve (b)
for a 2D lattice with linear size $n=8$. (c) Standard MPS representation for
1D chain. (d) Binary TTN Ansazt for 1D lattice. Yellow circles represent the
physical sites of the lattice, whereas the green circles indicate the tensors
that make up the TTN.

Beyond analytical approaches and exact solutions, that are typically limited
to a small class of Hamiltonians and microscopical interactions, numerical
methods play a key role in characterizing the low-energy properties of
strongly correlated quantum many-body systems, characterizing their phase
diagrams and the respective phase transitions [11, 12]. Despite their
impressive success and accuracy in many scenarios, the complete understanding
of several phenomena, e.g. the low-energy properties of Hubbard-like
Hamiltonians for interacting fermions or the behavior of spin systems with
frustrating interactions, is still debated, especially in two or three spatial
dimensions [13, 14, 15, 16, 17, 18, 19, 20], where current and near-future
experimental setups require benchmarking tools capable to simulate large
quantum many-body systems [21, 22, 23, 24].

In recent years, these challenging problems have inspired the development of
novel numerical algorithms based on Tensor Networks (TN), which do not suffer
from the well-known sign problem affecting Quantum Monte Carlo (QMC) methods
[25, 26, 27, 28]. These techniques exploit the area-law entanglement bounds
satisfied by a large class of physical systems and allow an efficient
representation of the low-energy sector by retaining in memory only the
relevant states that contribute to the equilibrium properties and the real-
time dynamics [29]. In particular, this is carried out by compressing the
exponentially large wavefunctions into a network of tensors interconnected
through auxiliary indices with bond dimension $m$. The main Anzätze for the
representation of quantum many-body states based on TNs include Matrix Product
States (MPS) for 1D systems [30, 31, 32], Projected Entangled Pair States
(PEPS) [33, 34, 35], Tree Tensor Networks (TTN) [36, 37, 38, 39] and
Multiscale Entanglement Renormalization Ansatz (MERA) [40, 41] which can be
defined in any dimension. However, while for one-dimensional systems MPS are
the established TN ansatz, the development of efficient TN algorithms for
higher-dimensional systems is still ongoing [42, 43, 44, 45, 46, 47, 48, 49].

In this framework, a possible and frequently used approach lies in the
application of consolidated one-dimensional numerical methods to the study of
higher-dimensional systems, by suitably mapping the high-dimensional problem
onto an equivalent one-dimensional one. For instance, the Density Matrix
Renormalization Group (DMRG) algorithm [50, 51], despite being defined for 1D
systems, has long been used for analysing two-dimensional quantum many-body
systems, providing important and reliable numerical results in a wide variety
of condensed-matter problems [52]. Two main strategies can be implemented: the
first one consists in considering the $n\times n^{\prime}$ 2D lattice as an
$n^{\prime}$-legs ladder, then grouping together all the $n^{\prime}$ sites of
each column in a single numerical site, and finally applying the one-
dimensional algorithms to the resulting chain [53]. The price to pay in this
case is the exponential growth with $n^{\prime}$ of the local basis dimension,
which usually limits the application of this strategy to quasi-two-dimensional
systems, in which the number of ladder legs is small with respect to the
longitudinal extension. The second strategy, on which we focus in this work,
lies in covering the two-dimensional lattice by a one-dimensional curve and
then applying one-dimensional algorithms to the resulting effective chain. The
drawback at this approach is that even if the original model contains only
nearest-neighbor interactions, the resulting model on the effective one-
dimensional chain shows long-range interactions, that have an influence on the
numerical approach efficiency [54].

In detail, the latter strategy induces a specific one-dimensional site-
ordering in the 2D lattice, generating a mapping between the 2D physical
system and a 1D chain. This procedure can be easily defined for any mapping

$\mathcal{M}:\,i\in[1,n]\times[1,n^{\prime}]\rightarrow\mu\in[1,nn^{\prime}]$
(1)

from a $n\times n^{\prime}$ lattice onto a chain with $n\,n^{\prime}$ sites,
through which it is possible to translate a two-dimensional model into a one-
dimensional one. Finally, one can apply 1D TN algorithms, such as DMRG, taking
into account the long-range interactions resulting from the mapping [55, 56,
57, 58, 59]. Long-range interactions eventually require a larger bond
dimension to properly describe the system. It follows that the performance of
the simulation depends on the capability of the mapping to preserve the
locality of the interactions from the original 2D Hamiltonian to the effective
new 1D one.

In this paper, we quantify the importance of the choice of the site-ordering
curve. In particular, we consider a square $n\times n$ lattice and compare the
precision we achieve by using the standard snake curve (see Fig. 1(a)) and the
Hilbert curve (see Fig. 1(b)) [60]: as we discuss in detail later, the
peculiar features of the Hilbert curve turn out to be very useful for
simulating systems of large sizes via TN methods, especially close to a
quantum phase transition where a divergence of the correlation length is
expected.

We present the numerical results for the 2D quantum Ising model in transverse
field obtained by using the Hilbert curve and the snake curve [61]. The
simulations are performed by using MPSs (see Fig. 1(c)) and binary TTNs (see
Fig. 1(d)). The latter are particularly interesting for a number of reasons
[39]: (i) with their hierarchical structure, they offer a strong connectivity
such that the distance between two lattice sites, in terms of links connecting
them, scales only logarithmically within the network, whereas for MPS the
distance is clearly linear. As a results, TTNs are particularly advantageous
to describe systems with long-range interactions with both open (OBC) and
periodic (PBC) boundary conditions; (ii) similarly to MPS, they provide an
efficient basis for the DMRG algorithm, which represent the state-of-the-art
technique for the numerical simulation of quantum many-body systems in 1D with
both local and long-range interactions; (iii) they have a loop-free structure
and, as a result, they show a low algorithmic complexity ($O(m^{4})$), that is
slightly larger than MPS ($O(m^{3})$), but significantly lower than intrinsic
2D TN algorithms such as PEPS or MERA, that suffer from very high
computational costs (at least $O(m^{10})$).

In order to test the performance obtained with the two mappings, we focus on
the ground state properties of the Ising model at zero temperature: we
numerically determine the ground state energy at different lattice sizes up to
$n=64$ with OBC. We find that the simulations performed with the Hilbert curve
mapping provide more accurate ground state energies, for both MPS and TTN
algorithms. We also find that the Hilbert curve shows a faster convergence of
the energy than the snake one as the bond dimension is increased. Finally, we
remarkably observe a larger improvement close to the critical point, stressing
the relevance of our analysis for the study of quantum phase transitions.
Analogous results are also provided for the PBC case with TTN and are shown in
the Appendix B.

The article is organized as follows. In Sec. I we introduce the reader to the
Hilbert and snake curves, briefly reviewing their space-filling and locality-
preserving properties. In Sec. II we describe the 2D quantum Ising model we
are going to focus on, while Sec. III is devoted to the presentation and
discussion of the numerical results we obtain, with the MPS and TTN
algorithms. In Sec. IV we draw our conclusions. Supplementary technical
details are given in the Appendix.

##  I The Hilbert Curve

The Hilbert curve is a fractal-like self-similar space-filling curve described
for the first time in 1891 by the mathematician David Hilbert [60]. Among
other interesting properties, it allows to create a mapping between a 1D and a
2D space by keeping and preserving the locality: this implies that two points
that are close to each other in the 1D space are close to each other also
after folding on the 2D space. The converse is not always strictly true, as
unavoidable when passing from two-dimensions to one-dimension. However, even
in this case, the curve shows a tendency to preserve the locality as much as
possible [62, 63, 64] which makes it a valuable tool for several applications
in computer science and bioinformatics [65, 66, 67]. The basic element of the
Hilbert curve is obtained by connecting the elements of a $2\times 2$ lattice
starting from the bottom-left (BL) to the bottom-right (BR) corner as shown in
Fig. 2(a). From this we can now construct the $n=4$ Hilbert curve via the
following procedure: draw the $n=2$ Hilbert curve into the top-left (TL) an
top-right (TR) quadrants, while rotating it clockwise and counterclockwise by
90 degrees in the BL and BR quadrants respectively. Then, join this different
curve replicas, ending up with the $n=4$ curve shown in Fig. 2(b).

(a)

(b)

(c)

(d)

Figure 2: (a) Hilbert curve for $n=2$, (b) for $n=4$, (c) for $n=8$. (d) Snake
curve for $n=8$. In (c) and (d) it is also shown the projection of a generic
site of the 2D lattice (within a given quadrant) and its four nearest
neighbours onto the one-dimensional chain.

The generalization to the level $n$ for covering the $n\times n$ lattice is
then straightforward: the $2n$ curve is obtained by drawing the $n$ curve in
the four main quadrants of the $2n\times 2n$ lattice and by applying the just
mentioned rules for the rotations in the BL and BR quadrants. In Fig. 2(c), by
way of example, we report the Hilbert curve for the $8\times 8$ lattice.

In the context of numerical simulations with TN algorithms, it would be
extremely useful to map a 2D lattice into a 1D chain in a way that preserves
locality at most, avoiding long-range interactions of sites that are very far
apart from each other. However, it is straightforward to prove that it is
impossible to map a $d$-dimensional lattice into a $d^{\prime}$-dimensional
one, with $d^{\prime}<d$, so that two neighboring sites in the $d$-lattice are
always close together in the $d^{\prime}$-lattice. Let us to consider, for
instance, nearest neighbour sites: each point of the $d$-dimensional square
lattice has $2d$ nearest neighbours, while only $2d^{\prime}$ nearest
neighbours in the $d^{\prime}$-dimensional one. Thus, after the mapping, there
will be at least $2(d-d^{\prime})$ sites placed at a distance larger than one
unit. For the specific case of a 2D square lattice and a 1D chain, this
implies that the optimal solution would be a mapping such that two out of four
nearest neighbours in the 2D lattice remain nearest neighbours in the 1D
chain. If we strictly follow this argument, the best choice would be the
frequently used snake curve, shown in Fig. 2(d). However, in this way, the
other two nearest neighbours of a generic site of the 2D lattice would be
rather far away from each other in the one-dimensional setup, with their
distance increasing up to $2n$. If our goal is to preserve the locality also
for these two other nearest neighbours, the curve that tends to globally
satisfy these constraints, outperforming other types of mapping, is exactly
the Hilbert curve. In order to get a general idea of this property, we report
in Fig. 2(c) the Hilbert mapping of a generic site of the 2D lattice (within a
given quadrant) and its four nearest neighbours: it is clearly visible that
all these points that are close on the two-dimensional plane remain fairly
close on the one-dimensional chain. This property, on the contrary, does not
hold true for the standard snake mapping, where, as expected, two out of four
nearest neighbours are almost always mapped at large distance, as shown in
Fig. 2(d). The only regions in which the Hilbert mapping does not preserve the
locality are the inner borders of the main quadrants. In this case, nearest
neighbor sites belonging to two adjacent quadrants turn out to be mapped quite
far away from each other.

##  II 2D Quantum Ising Model

In order to test the capabilities of the Hilbert curve mapping for our TN
algorithms, we consider the 2D quantum Ising model in presence of a uniform
transverse magnetic field along the $z$-direction [61]. For a $n\times n$
square lattice, the Hamiltonian reads

$H=J\sum_{\left<i,j\right>}\sigma_{i}^{x}\sigma_{j}^{x}+\lambda\sum_{i}\sigma_{i}^{z}\;,$
(2)

where $\left<i,j\right>$ refers to the nearest neighbour sites and
$\sigma_{i}^{\alpha}$ is the Pauli matrix in the $\alpha$-direction defined on
the site $i$. The coefficient $J$ represents the strength of the interaction
between nearest neighbour spins along the $x$-axis, whereas $\lambda$
determines the strength of the external magnetic field along the $z$-axis. In
the following, we consider the antiferromagnetic (AF) scenario for $J>0$, and,
without loss of generality, we fix the energy scale by setting $J=1$.

The two terms of the Hamiltonian in Eq. (2) do not commute since
$\left[\sigma_{i}^{x},\sigma_{i}^{z}\right]\neq 0$. Thus, there exist quantum
fluctuations in the system that can be tuned through the transverse field
term. For $\lambda=0$, the ground state shows a complete antiferromagnetic
order along the $x$-axis, producing a spontaneous staggered magnetization
$M_{s}=\frac{1}{n^{2}}\sum_{i}\left<GS|\zeta_{i}\sigma_{i}^{x}|GS\right>\neq
0$, in which $\zeta_{i}=\pm 1$ depending on the parity of the site $i$. In the
limit of $|\lambda|\to\infty$, on the contrary, each spin align with the
magnetic field along the $z$-axis, and, as a result, the antiferromagnetic
order along the $x$-direction vanishes, i.e., $M_{s}=0$. Between these two
regimes, the system is expected to undergo a quantum phase transition at
finite $|\lambda|$. In particular, we expect a second-order phase transition
at $|\lambda_{c}|\approx 3.1$ [68, 69].

Figure 3:  Staggered magnetization computed with TTN at OBC for different $n$
as a function of $\lambda$ and $m=25$. We observe that the magnetization drops
to zero close to the expected critical point. The gray line indicates the
expected position of the critical point.

We apply the MPS and the TTN algorithms to numerically compute the ground
state of the Hamiltonian in Eq. (2) at zero temperature. We have preliminarily
checked that our simulations capture the phase transition point by computing
the magnetization with the TTN algorithm, as shown in Fig. 3 for $n=4,8,16,32$
at OBC: As expected, the magnetization drops from one to zero close to the
transition point, whose expected position is indicated by the dashed gray line
[68, 69]. If not differently specified, all the results shown in the following
are obtained at OBC.

##  III Numerical results

Hereafter we compare the results obtained by exploiting the Hilbert and the
snake mapping with both the MPS and TTN wave function representations. The
Hilbert and the snake curves introduced above define two mappings
$\mathcal{M}=\mathcal{H},\mathcal{S}$ through which, from Eq. (2), we obtain
two different long-range, one-dimensional Hamiltonians

$H_{\mathcal{M}}=\sum_{<\mu,\nu>_{\mathcal{M}}}\sigma^{x}_{\mu}\sigma_{\nu}^{x}+\lambda\sum_{\mu}\sigma_{\mu}^{z}\,,$
(3)

where $<,>_{\mathcal{M}}$ indicates the sites mapped from nearest-neighbor
coordinates in the 2D lattice according to the mapping $\mathcal{M}$.

Figure 4:  Ground state energy density for the 2D Ising Hamiltonian (2) of a
MPS on Hilbert and snake mappings, evaluated for $n=4,8,16,32$ as a function
of the bond dimensions $m$ with OBC. Note how the energy difference due to the
use of the Hilbert with respect to the snake one increases with the size. Here
$\lambda=2.9$ and $J=1$.

In order to estimate the performance of a simulation with a given mapping and
representation, we compute the ground state energy density
$E_{\mathcal{H},\mathcal{S}}$ and the difference $\Delta
E=E_{\mathcal{S}}-E_{\mathcal{H}}$, with both the MPS and TTN algorithms. We
focus our analysis on the case $\lambda=2.9$ because, being close to
criticality, it offers the possibility to test the two mappings in a case in
which numerical simulations are expected to be more involving due to the
presence of large correlation length. Moreover, we have numerically verified
that this is the point where $\Delta E$ is larger with both MPS and TTN
representation, making our analysis clearer (for details see Appendix A).

First, we show the results obtained by employing the MPS Ansatz for computing
the Hamiltonian ground state. We perform simulations for the two effective
Hamiltonians $H_{\mathcal{H}}$ and $H_{\mathcal{S}}$ for $n=4,\ 8,\ 16,\ 32$,
for $\lambda=2.9$, and measure the ground state energy density (i.e. the
energy per site) as a function of the bond dimension $m$ (see Fig. 4). In each
panel, for both mappings we observe that the ground state energy density
decreases with the increase of the bond dimension, converging to an asymptotic
energy value. We notice that a clear improvement is achieved when the Hilbert
curve is adopted, with the energy difference increasing as the system size
grows. Analogous results are observed in Fig. 5 which have been obtained for
the same parameter values using now the TTN Ansazt with OBC: also in this case
the ground state energy densities associated with the Hilbert curve mapping
are smaller than those obtained with the snake curve. This behavior emerges
also for PBC, as shown in Appendix B.

Figure 5:  Ground state energy density for the 2D Ising Hamiltonian (2)
computed with TTN for $n=4,8,16,32$ as a function of the bond dimensions $m$
for OBC. Here $\lambda=2.9$ and $J=1$. As for the MPS case the energy
difference due to the use of the Hilbert with respect to the snake one
increases with the size.

In order to compare the results obtained with the MPS and the TTN Ansatze, we
first plot the ground state energy densities for the Hilbert curve case at
$m=50$. As shown in Fig. 6(a), we find that the values returned by the TTN
simulations are slightly lower than those found by exploiting MPS. We now
consider the difference $\Delta E$ as a function of $n$, plotted in Fig. 6(b)
up to $n=64$. We find that the Hilbert curve improves the energy estimation at
the second decimal digit, thus outperforming the snake mapping at each $n$
both with MPS and TTN. Moreover, while at low $n$ the improvement achieved
with the Hilbert mapping is larger with TTN, MPS exhibits a larger $\Delta E$
at higher values of $n$.

(a)

(b)

Figure 6:  (a) Ground state energy density computed with MPS and TTN by
exploiting the Hilbert mapping as a function of $n$ ($\lambda=2.9$, $m=50$).
(b) Energy density difference $\Delta E$ computed with the MPS and the TTN
algorithms as a function of $n$.

In order to provide a qualitative explanation to our findings, we have
analyzed the long range interaction terms generated by the two mappings, by
computing the probability distributions for the distances of the interaction
terms of the Hamiltonians $H_{\mathcal{H}}$ and $H_{\mathcal{S}}$. In
particular, we have measured the distances $d_{\mathrm{MPS}}$ and
$d_{\mathrm{TTN}}$, defined as the number of links connecting two different
lattice sites within the MPS and TTN network geometries respectively.

|
---|---
|

(a)

(b)

(c)

(d)

Figure 7:  Distributions of the distances $d_{\mathrm{MPS}}$ (first row) and
$d_{\mathrm{TTN}}$ (second row) relative to long range interaction terms in
the Hamiltonians $H_{\mathcal{H}}$ and $H_{\mathcal{S}}$ computed for $n=16$.

In Fig. 7 we plot the distributions of the distances $d_{\mathrm{MPS}}$ and
$d_{\mathrm{TTN}}$ with $n=16$ for the snake (left column) and the Hilbert
(right column) mappings: one can qualitatively see that the probabilities at
large distances are larger for the snake mapping than for the Hilbert one.
Furthermore, this feature is particularly evident in the TTN case, due to the
logarithmic scaling of distances within the network. It is possible to
understand which regions of the lattice contribute to the energy improvement
by studying the expectation values of the local magnetization. We focus on the
TTN simulations, measuring the local expectation values of the magnetization
$\\{\langle\sigma^{z}_{\mu}\rangle\\}_{1\leq\mu\leq n^{2}}$ on both the ground
states of the two Hamiltonians $H_{\mathcal{H}}$ and $H_{\mathcal{S}}$ and
computing the difference
$\Delta\sigma^{z}_{\mu}=|\langle\sigma^{z}_{\mu}\rangle_{\mathcal{H}}-\langle\sigma^{z}_{\mu}\rangle_{\mathcal{S}}|$,
as shown in Fig. 8. There are regions in which the difference is larger while
in others the magnetization values are similar. The former are located in the
bulk of the lattice main quadrants (regions delimited by black lines in Fig.
2(c)) where the Hilbert curve preserves the locality of interactions better
than the snake one. Instead, the latter are the regions close to the quadrant
boundaries, where the Hilbert curve fails to preserve locality and the
magnetization results are not more precise than those obtained with the snake
one.

To explain, instead, why the improvement observed at large $n$ is larger for
MPS with respect to the TTNs, we have checked how the fraction of interaction
terms with largest distances changes when we move from $n=16$ to $n=32$. We
have focused on the case of the Hilbert curve, finding that this fraction
decreases faster in the MPS case than in the TTN one, as shown in Fig. 9. For
this reason, we expect the improvement achievable with the MPS and the Hilbert
curve is larger than the one achievable with TTNs, even though the TTN
simulations are overall more accurate than the MPS ones.

|
---|---

(a)

(b)

Figure 8:  Local magnetization difference $\Delta\sigma^{z}_{\mu}$ between the
Hilbert and the snake ground states for $n=8$ (left panel) and $n=16$ (right
panel) computed by using TTNs for $\lambda=2.9,\,m=50$. The regions in which
the expectation values are similar are those in which neither the Hilbert nor
the snake curve are able to properly preserve the locality of the
interactions. .

|
---|---
|

(a)

(b)

(c)

(d)

Figure 9:  Distributions of the distances $d_{\mathrm{MPS}}$ (first row) and
$d_{\mathrm{TTN}}$ (second row) relative to the Hilbert mapping relative to
long range interaction terms computed for $n=16$ (first column) and $n=32$
(second column).

##  IV Conclusions

In this work, by using the MPS and the TTN algorithms, we have examined two
different strategies to study two-dimensional quantum many-body systems. Both
of them rely on the fact that it is possible to map a $d-$dimensional lattice
onto a one-dimensional one, by exploiting a suitable space-filling curve. In
particular, we have carried out a comparison between two possible mappings,
one based on the Hilbert curve and one on the snake ordering, finding a clear
advantage in the use of the former, both in the MPS and in TTN approach. While
the locality properties of the Hilbert ordering are directly mapped in the MPS
chain, in the case of TTNs they are enhanced by the logarithmic scaling of
distances within the TTN structure. By recalling that the standard way to use
TTNs for two-dimensional systems relies directly on the construction of a
binary 2D TTN [48], it is worth noting that the binary 2D TTN and the Hilbert
curve mapping generate the same long range interactions within the
corresponding networks. Therefore, the former benefits from the mathematical
locality-preserving properties of the Hilbert curve. From this point of view,
our analysis justifies the capability of binary 2D TTNs on the basis of the
interaction range of the effective 1D model induced by the network structure.

In conclusion, our work indicates a systematic strategy to evaluate, given a
quantum many-body system defined on a two-dimensional lattice, the more
efficient site-ordering to perform simulations by using TNs, based on the
analysis of the effective, one-dimensional system generated by the mapping. In
particular our work suggests that 1D variational ground state searching
algorithms, such as DMRG, strongly benefits from the choice of the Hilbert
curve for treating 2D systems. Other possible outlooks include the application
of this strategy to the study of systems characterized by non-trivial lattice
structures and interactions, as for example in spin glasses [70] and quantum
chemistry [71].

##  V Acknowledgments

We kindly acknowledge fruitful discussion with Timo Felser and Pietro Silvi.
MPS calculations were performed using the TeNPy Library [72]. This work is
partially supported by the Italian PRIN2017 and Fondazione CARIPARO, the INFN
project QUANTUM, the QuantERA projects QTFLAG and QuantHEP, the DFG project
TWITTER, the Horizon 2020 research and innovation programme under grant
agreement No 817482 (Quantum Flagship - PASQuanS). We acknowledge
computational resources by the Cloud Veneto, CINECA, the BwUniCluster, ATOS
Bull, and by University of Padova Strategic Research Infrastructure Grant
2017: “CAPRI: Calcolo ad Alte Prestazioni per la Ricerca e l’Innovazione”. AA,
NDP, and VG acknowledge support by MIUR via PRIN 2017 (Progetto di Ricerca di
Interesse Nazionale): project QUSHIP (2017SRNBRK).

##  Appendix A Choice of $\lambda=2.9$

We have computed the energy difference $\Delta E$ for both the MPS and the TTN
representatiosn with $m=20,35$ for values of $\lambda$ ranging in the interval
$[0,4]$ which includes the critical point. Then we have chosen, as explained
in the main text, the value of $\lambda$ where we observe a peak in the
precision improvement. Note that a clear improvement emerges at low bond
dimension values, implying that the Hilbert mapping guarantees a faster
convergence of the ground state energy. The results are shown in Fig. 10: We
have computed the energy difference $\Delta E$ at $n=8$ with OBC and PBC by
using the TTN algorithm (first row) and at $n=16$ with OBC by using the MPS
ground state representation (second row). In both rows, the left and right
panels refer to the results obtained with $m=20,35$, respectively.

|
---|---
|

(a)

(b)

(c)

(d)

Figure 10:  Difference $\Delta E$ of the ground state energy densities
computed with bond dimension $m=20$ (left column) and $m=35$ (right column).
In the first row TTN results are shown for $n=8$ with OBC and PBC, while in
the second row OBC MPS results are shown for $n=16$.

##  Appendix B Simulations with PBC

Figure 11:  Ground state energy density computed with TTN for $n=4,8,16$ as a
function of the bond dimensions $m$ for PBC. Note how the energy difference
due to the use of the Hilbert with respect to the snake one increases with the
size, as found in the OBC case.

Here we show the results obtained by employing the TTN Ansatz for the
computation of the ground state with the Hilbert and the snake mappings. As
anticipated in the main text, we observe an improvement when the Hilbert
mapping is adopted, as shown in Fig. 11.

## References

  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Bernevig and Hughes [2013] B. A. Bernevig and T. L. Hughes,  _Topological Insulators and Topological Superconductors_, stu - student edition ed. (Princeton University Press, 2013).
  * Girvin and Yang [2019] S. M. Girvin and K. Yang,  _Modern Condensed Matter Physics_ (Cambridge University Press, 2019).
  * Hansson  _et al._ [2017] T. H. Hansson, M. Hermanns, S. H. Simon, and S. F. Viefers, Quantum hall physics: Hierarchies and conformal field theory techniques, Rev. Mod. Phys. 89, 025005 (2017).
  * Sacépé  _et al._ [2020] B. Sacépé, M. Feigel’man, and T. M. Klapwijk, Quantum breakdown of superconductivity in low-dimensional materials, Nature Physics 16, 734 (2020).
  * Rachel [2018] S. Rachel, Interacting topological insulators: a review, Reports on Progress in Physics 81, 116501 (2018).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Reports on Progress in Physics 80, 016502 (2016).
  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, D. Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, The quantum technologies roadmap: a european community view, New Journal of Physics 20, 080201 (2018).
  * Alexeev  _et al._ [2020] Y. Alexeev, D. Bacon, K. R. Brown, R. Calderbank, L. D. Carr, F. T. Chong, B. DeMarco, D. Englund, E. Farhi, B. Fefferman, A. V. Gorshkov, A. Houck, J. Kim, S. Kimmel, M. Lange, S. Lloyd, M. D. Lukin, D. Maslov, P. Maunz, C. Monroe, J. Preskill, M. Roetteler, M. Savage, and J. Thompson, Quantum computer systems for scientific discovery, arXiv e-prints (2020), arXiv:1912.07577 [quant-ph] .
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * Scalapino [1998] D. Scalapino,  _Numerical Methods for Lattice Quantum Many-Body Problems_ (Addison-Wesley Longman, Incorporated, 1998).
  * Fehske  _et al._ [2008] H. Fehske, R. Schneider, and A. Weiße,  _Computational Many-Particle Physics_, Vol. 739 (2008).
  * Diep [2013] H. T. Diep,  _Frustrated Spin Systems_, 2nd ed. (World Scientific, 2013).
  * LeBlanc  _et al._ [2015] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Solutions of the two-dimensional hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5, 041041 (2015).
  * Qin  _et al._ [2020] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Absence of superconductivity in the pure two-dimensional hubbard model, Phys. Rev. X 10, 031016 (2020).
  * Matsuda  _et al._ [2013] Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Magnetization of ${\mathrm{srcu}}_{2}({\mathrm{bo}}_{3}{)}_{2}$ in ultrahigh magnetic fields up to 118 t, Phys. Rev. Lett. 111, 137204 (2013).
  * Stapmanns  _et al._ [2018] J. Stapmanns, P. Corboz, F. Mila, A. Honecker, B. Normand, and S. Wessel, Thermal critical points and quantum critical end point in the frustrated bilayer heisenberg antiferromagnet, Physical Review Letters 121, 10.1103/physrevlett.121.127201 (2018).
  * Maksimov  _et al._ [2019] P. Maksimov, Z. Zhu, S. R. White, and A. Chernyshev, Anisotropic-exchange magnets on a triangular lattice: Spin waves, accidental degeneracies, and dual spin liquids, Physical Review X 9, 10.1103/physrevx.9.021017 (2019).
  * Czarnik  _et al._ [2020] P. Czarnik, M. M. Rams, P. Corboz, and J. Dziarmaga, Tensor network study of the $m=1/2$ magnetization plateau in the shastry-sutherland model at finite temperature, (2020), arXiv:2012.05265 [cond-mat.str-el] .
  * Wietek  _et al._ [2020] A. Wietek, Y.-Y. He, S. R. White, A. Georges, and E. M. Stoudenmire, Stripes, antiferromagnetism, and the pseudogap in the doped hubbard model at finite temperature, (2020), arXiv:2009.10736 [cond-mat.str-el] .
  * Scholl  _et al._ [2020] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. LÃ€uchli, and A. Browaeys, Programmable quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms (2020), arXiv:2012.12268 [quant-ph] .
  * Ebadi  _et al._ [2020] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletic, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator (2020), arXiv:2012.12281 [quant-ph] .
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nat. Phys. 8, 277 (2012).
  * Henriet  _et al._ [2020] L. Henriet, L. Beguin, A. Signoles, T. Lahaye, A. Browaeys, G.-O. Reymond, and C. Jurczak, Quantum computing with neutral atoms, Quantum 4, 327 (2020).
  * Verstraete  _et al._ [2008] F. Verstraete, V. Murg, and J. Cirac, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Advances in Physics 57, 143 (2008), https://doi.org/10.1080/14789940801912366 .
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Unifying time evolution and optimization with matrix product states, Phys. Rev. B 94, 165116 (2016).
  * Montangero [2018] S. Montangero,  _Introduction to Tensor Network Methods_ (Springer International Publishing, 2018).
  * Orús [2019] R. Orús, Tensor networks for complex quantum systems, Nature Reviews Physics 1, 538 (2019), arXiv:1812.04011 [cond-mat.str-el] .
  * Eisert  _et al._ [2010] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Reviews of Modern Physics 82, 10.1103/revmodphys.82.277 (2010).
  * Fannes  _et al._ [1992] M. Fannes, B. Nachtergaele, and R. F. Werner, Finitely correlated states on quantum spin chains, Communications in Mathematical Physics 144, 443 (1992).
  * Klumper  _et al._ [1991] A. Klumper, A. Schadschneider, and J. Zittartz, Equivalence and solution of anisotropic spin-1 models and generalized t-J fermion models in one dimension, Journal of Physics A Mathematical General 24, L955 (1991).
  * Klümper  _et al._ [1993] A. Klümper, A. Schadschneider, and J. Zittartz, Matrix product ground states for one-dimensional spin-1 quantum antiferromagnets, EPL (Europhysics Letters) 24, 293 (1993).
  * Verstraete  _et al._ [2006] F. Verstraete, M. M. Wolf, D. Perez-Garcia, and J. I. Cirac, Criticality, the area law, and the computational power of projected entangled pair states, Phys. Rev. Lett. 96, 220601 (2006).
  * Verstraete and Cirac [2004] F. Verstraete and J. I. Cirac, Renormalization algorithms for Quantum-Many Body Systems in two and higher dimensions, (2004), arXiv:cond-mat/0407066 [cond-mat.str-el] .
  * Tepaske and Luitz [2020] M. Tepaske and D. J. Luitz, Three-dimensional isometric tensor networks, (2020), arXiv:2005.13592 [cond-mat.str-el] .
  * Shi  _et al._ [2006] Y.-Y. Shi, L.-M. Duan, and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Phys. Rev. A 74, 022320 (2006).
  * Gerster  _et al._ [2014] M. Gerster, P. Silvi, M. Rizzi, R. Fazio, T. Calarco, and S. Montangero, Unconstrained tree tensor network: An adaptive gauge picture for enhanced performance, Phys. Rev. B 90, 125154 (2014).
  * Silvi  _et al._ [2010] P. Silvi, V. Giovannetti, S. Montangero, M. Rizzi, J. I. Cirac, and R. Fazio, Homogeneous binary trees as ground states of quantum critical hamiltonians, Phys. Rev. A 81, 062335 (2010).
  * Silvi  _et al._ [2019] P. Silvi, F. Tschirsich, M. Gerster, J. Junemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes , 8 (2019).
  * Vidal [2007] G. Vidal, Entanglement renormalization, Phys. Rev. Lett. 99, 220405 (2007).
  * Evenbly and Vidal [2009] G. Evenbly and G. Vidal, Entanglement renormalization in two spatial dimensions, Phys. Rev. Lett. 102, 180406 (2009).
  * Phien  _et al._ [2015] H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing, Phys. Rev. B 92, 035142 (2015).
  * Vanderstraeten  _et al._ [2016] L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Gradient methods for variational optimization of projected entangled-pair states, Phys. Rev. B 94, 155123 (2016).
  * Fishman  _et al._ [2018] M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for contracting infinite two-dimensional tensor networks, Phys. Rev. B 98, 235148 (2018).
  * Xie  _et al._ [2017] Z. Y. Xie, H. J. Liao, R. Z. Huang, H. D. Xie, J. Chen, Z. Y. Liu, and T. Xiang, Optimized contraction scheme for tensor-network states, Phys. Rev. B 96, 045128 (2017).
  * Corboz  _et al._ [2018] P. Corboz, P. Czarnik, G. Kapteijns, and L. Tagliacozzo, Finite correlation length scaling with infinite projected entangled-pair states, Phys. Rev. X 8, 031031 (2018).
  * Kshetrimayum  _et al._ [2019] A. Kshetrimayum, M. Rizzi, J. Eisert, and R. Orús, Tensor network annealing algorithm for two-dimensional thermal states, Phys. Rev. Lett. 122, 070502 (2019).
  * Felser  _et al._ [2020] T. Felser, S. Notarnicola, and S. Montangero, Efficient tensor network ansatz for high-dimensional quantum many-body problems, (2020), arXiv:2011.08200 [quant-ph] .
  * Magnifico  _et al._ [2020] G. Magnifico, T. Felser, P. Silvi, and S. Montangero, Lattice quantum electrodynamics in (3+1)-dimensions at finite density with tensor networks, (2020), arXiv:2011.10658 [hep-lat] .
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * Schollwöck [2005] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Stoudenmire and White [2012] E. Stoudenmire and S. R. White, Studying two-dimensional systems with the density matrix renormalization group, Annual Review of Condensed Matter Physics 3, 111 (2012), https://doi.org/10.1146/annurev-conmatphys-020911-125018 .
  * Tschirsich  _et al._ [2019] F. Tschirsich, S. Montangero, and M. Dalmonte, Phase Diagram and Conformal String Excitations of Square Ice using Gauge Invariant Matrix Product States, SciPost Phys. 6, 28 (2019).
  * Yamada  _et al._ [2011] S. Yamada, T. Imamura, and M. Machida, Parallelization design on multi-core platforms in density matrix renormalization group toward 2-d quantum strongly-correlated systems, in  _SC ’11: Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis_ (2011) pp. 1–10.
  * Liang and Pang [1994] S. Liang and H. Pang, Approximate diagonalization using the density matrix renormalization-group method: A two-dimensional-systems perspective, Phys. Rev. B 49, 9214 (1994).
  * Xiang  _et al._ [2001] T. Xiang, J. Lou, and Z. Su, Two-dimensional algorithm of the density-matrix renormalization group, Phys. Rev. B 64, 104414 (2001).
  * Kantian  _et al._ [2019] A. Kantian, M. Dolfi, M. Troyer, and T. Giamarchi, Understanding repulsively mediated superconductivity of correlated electrons via massively parallel density matrix renormalization group, Phys. Rev. B 100, 075138 (2019).
  * Ghelli  _et al._ [2020] G. Ghelli, G. Magnifico, C. D. E. Boschi, and E. Ercolessi, Topological phases in two-legged heisenberg ladders with alternating interactions, Phys. Rev. B 101, 085124 (2020).
  * Ali [2021] M. Ali, On the ordering of sites in the density matrix renormalization group using quantum mutual information, (2021), arXiv:2103.01111 [quant-ph] .
  * Hilbert [1891] D. Hilbert, Ueber die stetige abbildung einer linie auf ein flÃ€chenstÃŒck, Mathematische Annalen 38, 459 (1891).
  * Suzuki  _et al._ [2013] S. Suzuki, J.-i. Inoue, and B. K. Chakrabarti,  _Quantum Ising Phases and Transitions in Transverse Ising Models_, Vol. 859 (2013).
  * Abel and Mark [1990] D. J. Abel and D. M. Mark, A comparative analysis of some two-dimensional orderings, International Journal of Geographical Information Systems 4, 21 (1990).
  * Jagadish [1990] H. V. Jagadish, Linear clustering of objects with multiple attributes, SIGMOD Rec. 19, 332 (1990).
  * Moon  _et al._ [2001] B. Moon, H. V. Jagadish, C. Faloutsos, and J. H. Saltz, Analysis of the clustering properties of the hilbert space-filling curve, IEEE Transactions on Knowledge and Data Engineering 13, 124 (2001).
  * Eavis and Cueva [2007] T. Eavis and D. Cueva, A hilbert space compression architecture for data warehouse environments, in  _Data Warehousing and Knowledge Discovery_ , edited by I. Y. Song, J. Eder, and T. M. Nguyen (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 1–12.
  * Lemire and Kaser [2011] D. Lemire and O. Kaser, Reordering columns for smaller indexes, Information Sciences 181, 2550 (2011).
  * Anders [2009] S. Anders, Visualization of genomic data with the Hilbert curve, Bioinformatics 25, 1231 (2009).
  * Friedman [1978] Z. Friedman, Ising model with a transverse field in two dimensions: Phase diagram and critical properties from a real-space renormalization group, Phys. Rev. B 17, 1429 (1978).
  * Li and Lei [2018] S.-H. Li and G.-P. Lei, Quantum phase transition in a two-dimensional quantum ising model: Tensor network states and ground-state fidelity, Journal of Physics: Conference Series 1087, 052011 (2018).
  * Mydosh [2015] J. A. Mydosh, Spin glasses: redux: an updated experimental/materials survey, Reports on Progress in Physics 78, 052501 (2015).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. D. Sawaya, S. Sim, L. Veis, and A. Aspuru-Guzik, Quantum chemistry in the age of quantum computing, Chemical Reviews 119, 10856 (2019).
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), SciPost Phys. Lect. Notes , 5 (2018), code available from https://github.com/tenpy/tenpy, arXiv:1805.00055 .


# Correlated Fractional Dirac Materials

Bitan Roy<EMAIL_ADDRESS>Department of Physics, Lehigh University,
Bethlehem, Pennsylvania, 18015, USA  Vladimir Juričić<EMAIL_ADDRESS>Departamento de Física, Universidad Técnica Federico Santa María, Casilla 110,
Valparaíso, Chile  Nordita, KTH Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, 106 91 Stockholm, Sweden

###### Abstract

Fractional Dirac materials (FDMs) feature a fractional energy-momentum
relation $E(\bm{k})\sim|\bm{k}|^{\alpha}$, where $\alpha\;(<1)$ is a real
noninteger number, in contrast to that in conventional Dirac materials with
$\alpha=1$. Here we analyze the effects of short- and long-range Coulomb
repulsions in two- and three-dimensional FDMs. Only a strong short-range
interaction causes nucleation of a correlated insulator that takes place
through a quantum critical point. The universality class of the associated
quantum phase transition is determined by the correlation length exponent
$\nu^{-1}=d-\alpha$ and dynamic scaling exponent $z=\alpha$, set by the band
curvature. On the other hand, the fractional dispersion is protected against
long-range interaction due to its nonanalytic structure. Rather a linear Dirac
dispersion gets generated under coarse graining and the associated Fermi
velocity increases logarithmically in the infrared regime, thereby yielding a
two-fluid system. Altogether, correlated FDMs unfolds a new landscape,
accommodating unconventional emergent many-body phenomena.

_Introduction_. Nodal Fermi liquids harbor a rich variety of electronic band
dispersions around a few isolated points in the Brillouin zone that are often
symmetry protected. They can display linear ($\alpha=1$), quadratic
($\alpha=2$) and cubic ($\alpha=3$) energy-momentum dispersion
$E(k_{j})\sim\pm|k_{j}|^{\alpha}$ along one or more than one component of the
spatial momentum ($\bm{k}$), for example. Here $+$ ($-$) corresponds to the
conduction (valence) band. Momentum is measured from the band touching point
at $\bm{k}=0$. As such, isotropic but linear and nonlinear band dispersions
can be observed in monolayer and multilayer graphene, respectively
graphene:RMP , and in three-dimensional Dirac and Weyl semimetals wsm:review .
Nodal Fermi liquids can also manifest mixed energy-momentum relations, where
integer $\alpha$ is direction dependent, as in multi-Weyl semimetals
multiweyl:1 ; multiweyl:2 ; multiweyl:3 ; multiweyl:4  and at the quantum
critical point (QCP) when the system resides at the brink of band insulation
SMInsQCP:1 ; SMInsQCP:2 ; SMInsQCP:3 . However, irrespective of microscopic
details, the symmetry class and the dimensionality of the system, $\alpha$ is
always an integer, hindering a large family of gapless electronic materials,
where $\alpha$ can take real fractional values. We name them _fractional Dirac
materials_ (FDMs), characterized by density of states (DOS)
$\varrho\sim|E|^{d/\alpha-1}$ and frequency ($\Omega$) dependent diagonal
optical conductivity (OC) $\sigma\sim\Omega^{(d-2)/\alpha}$ in the
noninteracting limit.

Figure 1: Constant energy contours in (a) $d=2$ for $\alpha=1$ (blue) and
$\alpha=2/3$ (red) (red) and (b) in $d=3$ for $\alpha=1$ (blue) and
$\alpha=2/3$ (yellow). We set $v_{\alpha}=1$ and $\bm{k}$ is measured in units
of $a^{-1}$, where $a$ is the lattice spacing [Eq. (1)]. Self-consistent
solutions of chiral symmetry breaking mass in (c) $d=2$ and (d) $d=3$ for
$\alpha=1$ (black), $0.8$ (blue), $0.6$ (red), $0.4$ (green) and $0.2$
(brown), as a function of the reduced distance from a QCP ($\delta$). For
$\alpha=1$ we recover conventional Dirac system with linearly dispersing
quasiparticles. Otherwise, the system describes an FDM of order $\alpha$. See
text for details.

Here we investigate the effects of short-range piece as well as of the long-
range tail of the repulsive Coulomb interactions in two-dimensional (2D) and
three-dimensional (3D) FDMs of order $\alpha<1$ for which the effective
Hamiltonian scales as $|k_{j}|^{\alpha}$ [Eq. (1) and Fig. 1]. Due to the
vanishing DOS in FDMs, sufficiently weak generic local or short-range
interactions are irrelevant perturbations. However, beyond a critical
threshold of interaction an FDM undergoes a continuous quantum phase
transition (QPT) through a QCP and becomes a correlated Dirac insulator. We
capture the university class of such QPT, which is distinct from its
counterpart in conventional Dirac systems ($\alpha=1$), within the framework
of a Gross-Neveu model gross-neveu . The associated critical exponents in turn
also govern the scaling of the spectral gap in the ordered phase with the
interaction strength [Fig. 1]. We show that the QPT in 2D (3D) FDMs is non-
Gaussian (mean-field or Gaussian) in nature.

In the presence of long-range Coulomb interaction, the electronic dispersion
in 2D and 3D FDMs, on the other hand, remains invariant due to its nonanalytic
structure [Eq. (1)]. Rather under coarse graining, as the system approaches
deep infrared regime, a linear Dirac dispersion gets generated through the
quantum many-body corrections. The effective Fermi velocity of such emergent
conventional Dirac quasiparticles then increases logarithmically at lower
energies. Resultant system thus describes a two-component quantum fluid at low
energies, composed of (a) effectively noninteracting (protected by
nonanalyticity of the dispersion) fractional and (b) marginally sharp (due to
logarithmically increasing Fermi velocity) conventional Dirac fermions.

Altogether, correlated fractional Dirac liquids unfold a previously unexplored
territory of rich emergent quantum phenomena, which can be studied using
quantum Monte Carlo simulation of their lattice realizations in terms of
infinitely long-ranged power-law hopping along the principal directions
syzranov:longrange  and possibly in electronic fractal lattices, recently
realized on designer quantum materials kempkes-fractal-exp ; shang-fractal-exp
. Fractal lattices are promising platforms to host FDMs, where the
characteristic irrational fractal dimension $d_{\rm frac}$ may give rise to
noninteger $\alpha$.

_Model_. The effective Hamiltonian for a $d$-dimensional FDM of order $\alpha$
takes the form

$~{}H_{\rm
FD}(\bm{k})=\sum^{d}_{j=1}v_{\alpha}|k_{j}|^{\alpha}\;\text{sign}(k_{j})\;\Gamma_{j},$
(1)

where $v_{\alpha}$ bears the dimension of energy $\times\;a^{\alpha}$, such
that $H_{\rm FD}(\bm{k})$ has the dimension of energy. We set the lattice
spacing $a=1$. Mutually anticommuting $D$-dimensional Hermitian $\Gamma$
matrices, operating on the orbital or sublattice degrees of freedom, satisfy
the algebra $\\{\Gamma_{j},\Gamma_{k}\\}=2\delta_{jk}\mathrm{I}_{D}$, where
$j,k=1,\cdots,D$ and $\mathrm{I}_{D}$ is a $D$-dimensional identity matrix.
For now we keep $D$ to be arbitrary. The energy spectra of $H_{\rm FD}$ are
$\pm E_{\alpha}(\bm{k})$, where

$E_{\alpha}(\bm{k})=v_{\alpha}\bigg{[}\;\sum^{d}_{j=1}|k_{j}|^{2\alpha}\;\bigg{]}^{\frac{1}{2}}\equiv
v_{\alpha}|k|^{\alpha}\bigg{[}\;\sum^{d}_{j=1}|\hat{\Omega}_{j}|^{2\alpha}\;\bigg{]}^{\frac{1}{2}}$
(2)

and $\hat{\Omega}_{j}$ are the components of a $d$-dimensional spherical unit
vector. Noninteracting $d$-dimensional FDMs are respectively characterized by
the following power-law scaling of DOS and diagonal OC supplementary

$\displaystyle~{}\varrho=D_{d}(\alpha)\;|E|^{\frac{d}{\alpha}-1}\>\;\text{and}\>\;\sigma=\frac{e^{2}}{h}\frac{\pi
D}{16}C_{d}(\alpha)\left(\frac{\Omega}{v_{\alpha}}\right)^{\frac{d-2}{\alpha}},$
(3)

The scaling of $D_{d}(\alpha)$ and $C_{d}(\alpha)$ are shown in Fig. 2. Thus
specific heat in FDMs scales as $C_{v}\sim T^{d/\alpha}$ with temperature
($T$).

The imaginary time ($\tau$) Euclidean action associated with $H_{\rm
FD}(\bm{k})$ reads

$S_{0}=\int d\tau\int
d^{d}\bm{x}\left\\{\Psi^{\dagger}\;\left[\partial_{\tau}+H_{\rm
FD}(\bm{k}\rightarrow-i\bm{\nabla})\right]\;\Psi\right\\},$ (4)

where we set $\hbar=1$, and $\Psi$ and $\Psi^{\dagger}$ are $D$-dimensional
independent Grassmann variables. Under coarse graining $\tau\to e^{z\ell}\tau$
and $\bm{x}\to e^{\ell}\bm{x}$, where $\ell$ is the logarithm of the
renormalization group (RG) scale and $z$ is the dynamic scaling exponent,
measuring the relative scaling between the energy and momentum, with
$z=\alpha$ in noninteracting FDMs. Then the scale invariance of the free
action $S_{0}$ mandates the following scaling dimensions
$[\Psi]=[\Psi^{\dagger}]=d/2$ and $[v_{\alpha}]=z-\alpha$. Hence, in a
noninteracting system $v_{\alpha}$ is scale invariant or marginal as
$\alpha=z$ therein. Next we seek to scrutinize the effects of interactions
among fractional Dirac quasiparticles. In what follows, we discuss the
imprints of short-range and long-range Coulomb interactions on FDMs
separately.

Figure 2: Scaling of two universal functions $D_{d}$ [(a) and (b)] and $C_{d}$
[(c) and (d)], respectively governing the energy and frenqency dependances of
DOS and OC in two ($d=2$) and three ($d=3$) dimensions. See Eq. (3).

_Short-range interaction_. Short-range interactions among electronic
quasiparticles in nodal Fermi liquids typically give birth to various
spontaneously broken symmetry phases. Among them Dirac masses are most
prominent at low temperatures. In the ordered phase, they give rise to
isotropic gapped quasiparticle spectra, yielding electrical or thermal
insulators, and lead to a maximal gain of the condensation energy. A Dirac
mass is represented by the fermion bilinear $\Phi=\Psi^{\dagger}M\Psi$. The
$D$-dimensional Hermitian matrix $M$ satisfies $\\{M,\Gamma_{j}\\}=0$ for
$j=1,\cdots,d$ and $M^{2}=\mathrm{I}_{D}$. In the ordered phase
$\langle\Psi^{\dagger}M\Psi\rangle\neq 0$. In the spirit of the Gross-Neveu
formalism such a Dirac mass can be favored by a local or momentum independent
four-fermion interaction for which the Euclidean action is gross-neveu

$S^{\rm SR}_{\rm int}=\int d\tau\int
d^{d}\bm{x}\;g_{{}_{m}}\left(\Psi^{\dagger}M\Psi\right)^{2},$ (5)

where $g_{{}_{m}}$ is the corresponding coupling constant. For simplicity,
here we consider a single-component microscopic Isinglike symmetry breaking
scalar Dirac mass that yields an insulator in the ordered phase. The following
discussion can be generalized to vectorlike Dirac masses, which, on the other
hand, break a continuous symmetry following the spirit of the Nambu-Jona-
Lasinio NJL  model as well as for superconducting Dirac masses roy-juricic:SC
. Leaving these cases for future investigations, here we focus on the
nucleation of an Isinglike Dirac mass.

The scale invariance of $S^{\rm SR}_{\rm int}$ leads to
$[g_{{}_{m}}]=\alpha-d$. Therefore, when $\alpha<d$, yielding a vanishing DOS
in FDMs, a sufficiently weak generic local quartic interaction is an
_irrelevant_ parameter. Consequently, nucleation of any Dirac mass takes place
beyond a critical threshold of interaction via a QPT occurring through a QCP
located at $g_{{}_{m}}=g^{\ast}_{{}_{m}}$ (say). Here we capture such an
emergent quantum critical phenomenon in a correlated FDM from a perturbative
RG analysis controlled by a small parameter $\epsilon=\alpha-d$. Upon
accounting for the leading order or one-loop quantum corrections we arrive at
the RG flow equation for the dimensionless coupling constant
$\lambda_{m}=g_{{}_{m}}\Lambda^{d-\alpha}S_{d}(\alpha)/v_{\alpha}$, explicitly
given by supplementary

$~{}\beta_{\lambda_{m}}=\frac{d\lambda_{m}}{d\ell}=-\epsilon\lambda_{m}+(D-2)\lambda^{2}_{m}$
(6)

after integrating out fast Fourier modes with Mastubara frequencies
$-\infty<\omega<\infty$ and residing within a thin Wilsonian momentum shell
$\Lambda e^{-\ell}<|\bm{k}|<\Lambda$. Here $\Lambda$ is the ultraviolet
momentum cutoff up to which the FDMs show fractional dispersion with power
$\alpha$ [Eq. (1)] and

$~{}S_{d}(\alpha)=\int\frac{d\hat{\Omega}}{(2\pi)^{d}}\bigg{[}\sum^{d}_{j=1}|\hat{\Omega}_{j}|^{2\alpha}\bigg{]}^{-\frac{1}{2}}\equiv\int\frac{d\hat{\Omega}}{(2\pi)^{d}}\frac{1}{[f(\hat{\Omega})]^{1/2}}.$
(7)

The scaling of $S_{d}(\alpha)$ in $d=2$ and $d=3$ are shown in Fig. 3. Notice
that when the $\Gamma$ matrices are Pauli matrices ($D=2$), the perturbative
corrections $\sim\lambda^{2}_{m}$ in Eq. (6) vanish. Similar conclusion holds
for conventional Dirac fermions up to the two-loop order gracey:twocomponent ,
which, however, possibly undergo a continuous QPT at finite coupling rosa-
vitale-witerich . Whether such a conclusion holds for FDMs remains to be
investigated.

Thus, let us consider $D=4$, the minimal dimensionality of the $\Gamma$
matrices for which the quantum corrections in the RG flow equation of
$\lambda_{m}$ [Eq. (6)] is nontrivial. Then 2D fractional and conventional
Dirac systems enjoy a continuous global SU(2) $\otimes$ U(1) chiral symmetry,
respectively generated by $\\{\Gamma_{34},\Gamma_{45},\Gamma_{35}\\}$ and
$\Gamma_{12}$, where $\Gamma_{jk}=[\Gamma_{j},\Gamma_{k}]/(2i)$. The Dirac
mass breaks the continuous SU(2) chiral symmetry, unless $M=\Gamma_{12}$. This
is so because the maximal number of mutually anticommuting four-component
Hermitian matrices is five. By contrast, in $d=3$ the Dirac mass breaks a
continuous U(1) chiral symmetry, generated by $\Gamma_{45}$.

Figure 3: Scaling of $S_{d}$ [(a) and (b)] and $F_{d}$ [(c) and (d)]
respectively appearing in the definition of dimensionless Gross-Neveu
interaction [Eq. (6)] and the RG flow equation of Fermi velocity
($v_{{}_{1}}$) of emergent linearly dispersing Dirac quasiparticles due to
long-range Coulomb repulsion [Eq. (12)] in two ($d=2$) and three ($d=3$)
dimensions.

The QCP describing the QPT between a nodal fractional Dirac liquid and a Dirac
insulator is located at $\lambda_{m}=\lambda^{\ast}_{m}=\epsilon/(D-2)$. As
the fermionic self-energy correction due to local interactions vanishes to the
leading order in the $\epsilon$ expansion supplementary , this QCP is
characterized by the correlation length exponent ($\nu$), given by

$\nu^{-1}=\left.\frac{d\beta_{\lambda_{m}}}{d\lambda_{m}}\right|_{\lambda_{m}=\lambda^{\ast}_{m}}=\epsilon=d-\alpha,\>\>\text{and}\>\>z=\alpha.$
(8)

These two exponents determine the universality class of the fractional Dirac
liquid to insulator QPT. In addition, they determine the scaling of the mass
gap, which can be demonstrated by solving the self-consistent gap equation.

The self-consistent gap equation is obtained by performing a Hubbard-
Stratonovich decomposition of the four-fermion term via a bosonic field $\Phi$
and subsequently integrating out the fermionic fields, yielding supplementary

$\frac{1}{g_{{}_{m}}}=\int^{\prime}\frac{d^{d}\bm{k}}{(2\pi)^{d}}\frac{1}{\left[E^{2}_{\alpha}(\bm{k})+\Phi^{2}\right]^{1/2}}\equiv
F(\Phi).$ (9)

The momentum integral is restricted up to the ultraviolet cutoff $\Lambda$
(denoted by the ‘prime’ symbol). The right hand side of this gap equation
scales as $\Lambda^{d-\alpha}$ and it is ultraviolet divergent for $d>\alpha$.
Such an ultraviolet divergence can be regularized by defining a critical
coupling for the ordering $g^{\ast}_{{}_{m}}=[F(0)]^{-1}$, in terms of which
we arrive at the regularized gap equation

$~{}\delta=\int^{1}_{0}x^{d-1}\left[\frac{1}{x^{\alpha}}-\int\frac{d\hat{\Omega}}{(2\pi)^{d}}\frac{[S_{d}(\alpha)]^{-1}}{[x^{2\alpha}f(\hat{\Omega})+m^{2}]^{1/2}}\right]dx,$
(10)

where $x=k/\Lambda$ and $m=\Phi/(v_{\alpha}\Lambda^{\alpha})$ are
dimensionless, and the reduced distance from the QCP is
$\delta=(\lambda_{m}-\lambda^{\ast}_{m})/(\lambda_{m}\lambda^{\ast}_{m})$. We
numerically solve this gap equation in $d=2$ and $d=3$. The results are shown
in Fig. 1.

Notice that nontrivial solution for the mass gap ($m$), and concomitantly for
the order parameter, exists only when $\delta>0$ or
$\lambda_{m}>\lambda^{\ast}_{m}$, i.e., when interaction strength
($\lambda_{m}$) is above a critical one ($\lambda^{\ast}_{m}$). The scaling of
$m$ with $\delta$ in FDMs is distinct from its counterparts in conventional
Dirac systems, where the mass gap shows linear (as $\nu=z=1$) and square-root
(as $2\nu=z=1$) scaling with $\delta$ in $d=2$ and $d=3$, respectively, since
$m\sim\delta^{\nu z}$. In 3D Dirac system the scaling of $m$ with $\delta$
displays a logarithmic correction due to the violation of the hyperscaling
hypothesis, as the system then lives at the upper critical dimension $d_{\rm
up}=3$ roydassarma . The upper critical dimension for FDMs is $d_{\rm
up}=2+\alpha$, where $\nu=1/2$. Therefore, 2D FDMs always remain below the
upper critical dimension and the QPT is non-Gaussian in nature. By contrast, a
3D FDM always lives above the upper critical dimension, and the QPT is
Gaussian or mean-field in nature.

_Long-range Coulomb interaction_. The instantaneous long-range Coulomb
interaction is known to renormalize the Fermi velocity of linearly dispersing
Dirac quasiparticles in both $d=2$ gonzalez-NPB  and $d=3$ goswami-PRL ;
hosur-PRL ; roy-juricic:longrange3D . This interaction is represented by a
scalar gauge field $a_{0}(\tau,\bm{x})$ minimally coupled to the density of
the quasiparticles. The corresponding propagator exhibits a dimensionality-
dependent scaling with momentum $\sim|\bm{k}|^{d-1}$, ensuring the
characteristic $1/r$ behavior of the density-density Coulomb interaction in
real space in any $d$. The form of the Coulomb propagator, in turn, dictates
that the electric charge ($e_{d}$) can (cannot) receive perturbative
corrections in $d=3$ ($d=2$). The Coulomb part of the Euclidean action
therefore reads as

$S_{\rm C}=\int d\tau\int
d^{d}\bm{x}\left(-ie_{d}a_{0}\Psi^{\dagger}\Psi+a_{0}\frac{1}{2|{\bm{\nabla}}|^{d-1}}a_{0}\right).$
(11)

To analyze the effect of the long-range Coulomb interaction on fractional
Dirac fermions, we compute the one-loop self-energy diagram. We find that the
parameter $v_{\alpha}$ appearing in Eq. (1) does not renormalize when
$\alpha\neq 1$, as a consequence of the nonanalytic structure of the
dispersion of fractional Dirac excitations. However, the long-range Coulomb
interaction generates linearly dispersing quasiparticles, for which the
Hamiltonian is given by Eq. (1) with $\alpha=1$. The RG flow equation for the
Fermi velocity ($v_{{}_{1}}$) of emergent linearly dispersing quasiparticles
in $d=2$ and $d=3$ reads supplementary

$\frac{dv_{{}_{1}}}{d\ell}=\frac{1}{C_{d}}\,F_{d}(\alpha)\,\alpha_{{}_{\rm
FS}}\,v_{{}_{1}},$ (12)

where $C_{2}=8\pi$ and $C_{3}=6\pi^{2}$, and the function (Fig. 3)

$F_{d}(\alpha)=\alpha\int\frac{d\hat{\Omega}}{N_{d}}\;\frac{|\hat{\Omega}_{i}|^{\alpha-1}}{[f(\hat{\Omega})]^{1/2}}\;\bigg{(}1-\frac{|\hat{\Omega}_{i}|^{2\alpha}}{f(\hat{\Omega})}\bigg{)},$
(13)

with $N_{2}=\pi$ and $N_{3}=8\pi/3$. The associated effective fine structure
constant is $\alpha_{{}_{\rm FS}}={e_{d}^{2}}/{v_{{}_{1}}}$. The function
$f(\hat{\Omega})$ is defined in Eq. (7), and $F_{d}(\alpha)$ is independent of
the choice of the component ($i=1,\cdots,d$) of the $d$-dimensional unit
vector ($\hat{\Omega}_{i}$) at least in $d=2$ and $3$.

Notice that the Fermi velocity $v_{{}_{1}}$ grows logarithmically as the
system approaches the deep infrared regime under coarse graining, as
$F_{2}(\alpha)$ and $F_{3}(\alpha)$ are both positive definite. Furthermore,
as the power of the fractional dispersion increases for $0<\alpha\leq 1$, both
functions increase monotonically, as shown in Fig. 3. This is, indeed,
consistent with the scaling of the DOS with energy. Namely, as the DOS
increases with increasing $\alpha$, the flow of the Fermi velocity of
generated Dirac quasiparticles also increases. When the bare quasiparticles
are linearly dispersing, well-known results for the flow of the Fermi
velocities in $d=2$ and $d=3$ are readily recovered gonzalez-NPB ; goswami-PRL
; hosur-PRL ; roy-juricic:longrange3D . Thus, long range coulomb interaction
in FDMs gives birth to a two-component quantum fluid constituted by (a)
effectively noninteracting fractional and (b) marginal (due to logarithmically
increasing Fermi velocity) conventional Dirac quasiparticles.

Finally, we show that in $d=3$ the long-range tail of the Coulomb interaction
is screened by fractional Dirac excitations. To this end, we compute the
polarization (bubble) diagram, and find that the Coulomb charge in $d=3$ is
logarithmically decreasing under coarse graining, with the RG flow equation of
the form

$\frac{de^{2}_{d}}{d\ell}=-\alpha_{\rm
FS}^{(\alpha)}\,e^{2}_{d}\,\mathcal{F}(\alpha)\>\>\text{or}\>\>\frac{d\alpha_{\rm
FS}^{(\alpha)}}{d\ell}=-\left[\alpha_{\rm
FS}^{(\alpha)}\right]^{2}\,\mathcal{F}(\alpha).$ (14)

Here $\alpha_{\rm FS}^{(\alpha)}={e^{2}_{d}}/(v_{\alpha}\Lambda^{\alpha-1})$
is the effective dimensionless fine structure constant of the FDM and the
function

$\mathcal{F}(\alpha)=\frac{\alpha}{4}\big{[}(1-3\alpha)I(\alpha)+2\alpha
J(\alpha)\big{]},$ (15)

for $1/2<\alpha\leq 1$, with $I(\alpha)$ and $J(\alpha)$ defined in terms of
the components of the three-dimensional unit vector

$\displaystyle I_{\rho\sigma}(\alpha)$
$\displaystyle=\int\frac{d\hat{\Omega}}{(2\pi)^{3}}\;\frac{|\hat{\Omega}_{\rho}|^{2\alpha-2}|\hat{\Omega}_{\sigma}|^{2\alpha-2}\hat{\Omega}_{\rho}\hat{\Omega}_{\sigma}}{[f(\hat{\Omega})]^{5/2}}\equiv
I(\alpha)\delta_{\rho\sigma},$ $\displaystyle J_{\rho}(\alpha)$
$\displaystyle=\int\frac{d\hat{\Omega}}{(2\pi)^{3}}\;\frac{|\hat{\Omega}_{\rho}|^{2\alpha-2}}{[f(\hat{\Omega})]^{3/2}}\equiv
J(\alpha).$ (16)

The function $\mathcal{F}(\alpha)$ is strictly positive, implying that the
charge $e_{d}$ and fine structure constant $\alpha_{\rm FS}^{(\alpha)}$ of the
FDM decrease during the RG flow. Furthermore, the fine structure constant of
generated Dirac quasiparticles ($\alpha_{{}_{\rm FS}}$) also decreases but
even faster than in a conventional relativistic Dirac material due to
additional screening by microscopic fractional Dirac excitations.

_Summary and discussions_. Here we explore the emergent quantum critical
behavior of correlated fractional Dirac liquids in two and three dimensions in
the presence of short-range and long-range Coulomb interactions. Strong short-
range interactions can give rise to spontaneous symmetry breaking. By
contrast, the long-range component of Coulomb interaction gives birth to
linear Dirac dispersion through a spectral reconstruction, while leaving the
original fractional Dirac dispersion unaffected, thereby yielding a two-
component quantum fluid in the infrared regime. In future, we will develop the
quantum critical description for the spontaneous symmetry breaking in terms of
strongly coupled boson-fermion Mott-Yukawa theory zinnjustin:book , which will
allow us to capture hallmarks of emergent non-Fermi liquid near the associated
QCP and the role of retarded boson-fermion Yukawa interaction on possible
restoration of the Lorentz symmetry in correlated FDMs RJH:JHEP . Furthermore,
by incorporating retarded current-current interaction, accompanying the
instantaneous density-density Coulomb repulsion, we will develop quantum
electrodynamics for FDMs.

_Acknowledgments_. B.R. was supported by the Startup grant from Lehigh
University. V.J. acknowledges support of the Swedish Research Council (VR
2019-04735). Nordita is partially supported by Nordforsk.

## References

  * (1) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * (2) N. P. Armitage, E. J. Mele, E. J. and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * (3) G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Phys. Rev. Lett. 107, 186806 (2011).
  * (4) C. Fang, M. J. Gilbert, X. Dai, and B. A. Bernevig, Phys. Rev. Lett. 108, 266802 (2012).
  * (5) B.-J. Yang and N. Nagaosa, Nat. Commun. 5, 4898 (2014).
  * (6) B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, aaf5037 (2016).
  * (7) G. Montambaux, F. Piéchon, J.-N. Fuchs, and M. O. Goerbig, Phys. Rev. B 80, 153412 (2009).
  * (8) B. Roy, P. Goswami and V. Juričić, Phys. Rev. B 95, 201102(R) (2017).
  * (9) B. Roy and M. P. Foster, Phys. Rev. X 8, 011049 (2018).
  * (10) D. J. Gross and A. Neveu, Phys. Rev. D 10, 3235 (1974).
  * (11) M. Garttner, S. V. Syzranov, A. M. Rey, V. Gurarie, and L. Radzihovsky, M. Garttner, Phys. Rev. B 92, 041406(R) (2015).
  * (12) S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J.M. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. M. Smith, Nat. Phys. 15, 127 (2019).
  * (13) J. Shang, Y. Wang, M. Chen, J. Dai, X. Zhou, J. Kuttner, G. Hilt, X. Shao, J. M. Gottfried, and K. Wu, Nat. Chem. 7, 389 (2015).
  * (14) See Supplementary Materials at XXX-XXXX for the derivation of DOS and OC, and details of RG analysis in FDMs with short- and longe-range Coulomb repulsion.
  * (15) Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961).
  * (16) B. Roy and V. Juričić, Phys. Rev. B 90, 041413(R) (2014).
  * (17) J. A. Gracey, Int. J. Mod. Phys. A 9, 727 (1994).
  * (18) L. Rosa, P. Vitale, and C. Wetterich, Phys. Rev. Lett. 86, 958 (2001).
  * (19) B. Roy and S. Das Sarma, Phys. Rev. B 94, 115137 (2016).
  * (20) J. Gonzalez, F. Guinea, and M. A. H. Vozmediano, Nucl. Phys. B 424, 595 (1994).
  * (21) P. Goswami and S. Chakravarty, Phys. Rev. Lett. 107, 196803 (2011).
  * (22) P. Hosur, S. A. Parameswaran, and A. Vishwanath, Phys. Rev. Lett. 108, 046602 (2012).
  * (23) B. Roy and V. Juričić, Phys. Rev. B 96, 155117 (2017).
  * (24) J. Zinn-Justin, _Quantum Field Theory and Critical Phenomena_ (Oxford University Press, Oxford, U.K., 2002).
  * (25) B. Roy, V. Juričić, and I. F. Herbut, J. High Energy Phys. 2016, 18 (2016).