# Optimal free descriptions of many-body theories

Christopher J. Turner<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Leeds, Leeds LS2 9JT, United Kingdom  Konstantinos
Meichanetzidis<EMAIL_ADDRESS>School of Physics and Astronomy, University
of Leeds, Leeds LS2 9JT, United Kingdom  Zlatko Papić  School of Physics and
Astronomy, University of Leeds, Leeds LS2 9JT, United Kingdom  Jiannis K.
Pachos  School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT,
United Kingdom

###### Abstract

Interacting bosons or fermions give rise to some of the most fascinating
phases of matter, including high-temperature superconductivity, the fractional
quantum Hall effect, quantum spin liquids and Mott insulators. While these
systems are promising for technological applications, they also present
conceptual challenges as they require approaches beyond mean-field and
perturbation theory. Here we develop a general framework for identifying the
free theory that is closest to a given interacting model in terms of their
ground state correlations. Moreover, we quantify the distance between them
using the entanglement spectrum. When this interaction distance is small, the
optimal free theory provides an effective description of the low energy
physics of the interacting model. Our construction of the optimal free model
is non-perturbative in nature, thus it offers a new theoretical framework for
investigating strongly correlated systems.

Many-body physics of non-interacting systems reduces to the description of a
single particle, and thus is well understood. In contrast, the understanding
of interacting systems remains one of the major open problems of both
condensed matter and high energy physics. Complete analytical solutions are
possible in the so-called integrable systems Sutherland , which are not robust
to perturbations. More generally, one relies on mean-field approaches, density
functional theory, or perturbation theory to expand around known solvable
instances. Such methods can be employed when correlations are weak, or
interactions induce small corrections to the original state. Many interesting
phenomena, however, have non-perturbative origin, such as superconductivity or
the fractional quantum Hall effect. While important insights about such
systems have been obtained using variational ansätze FeynmanRoton ; BCS ;
Gutzwiller ; Laughlin , this approach requires non-trivial physical intuition
about the nature of the emerging free quasiparticles. A question of paramount
importance arises: is it possible to directly identify the free effective
model that is most similar to a given interacting one?

Here we introduce the concept of the interaction distance, $D_{\cal F}$, which
quantifies the effect of interactions on the ground state of a many-body
system. At the same time we identify the optimal free theory closest to the
given interacting model. Our approach employs quantum information inspired
techniques to study the correlations of a system witnessed in its entanglement
spectrum and to build a general diagnostic tool of interactions. Typically we
find $D_{\mathcal{F}}$ to be small when mean-field theory is applicable, while
non-trivial behaviour emerges near critical regions. Using the example of the
quantum Ising model, we demonstrate that the interaction distance can be
calculated efficiently. We envision that finding the optimal free description
of interacting systems can help formulate suitable variational ansätze in a
variety of areas, ranging from condensed matter to high energy physics.
Alternatively, it could be used to develop efficient numerical simulations of
interacting systems that scale favourably with the system size.

##  I Results

Overview:– The interaction distance $D_{\cal F}$, introduced below, quantifies
the distance between the reduced density matrix of a many-body state and the
closest free-particle reduced density matrix. We demonstrate that $D_{\cal F}$
can be calculated efficiently from the entanglement spectrum Haldane08 , and
when this distance is small ($D_{\cal F}\ll 1$) leads to the optimal free
model that best describes the low energy properties of the interacting system.
Near critical regions, $D_{\cal F}$ behaves non-trivially, and its finite-size
scaling can be related to the properties of the model under renormalisation
group flow. As an example, we apply our method to the 1D quantum Ising model
in the presence of transverse and longitudinal fields. We demonstrate that
this model has $D_{\cal F}\approx 0$ across the whole phase diagram and we
identify its optimal free description. Finally, we present a particular model
with non-zero interaction distance in the thermodynamic limit, thus
demonstrating its intrinsic interacting nature.

Interaction distance and optimal free model:– Consider an arbitrary many-body
system prepared in its ground state, $\left|\Psi\right\rangle$. For
simplicity, here we consider fermionic systems defined on a lattice, although
our approach can be generalised to other systems, as we discuss below.
Partitioning the system in two regions, $A$ and its complement $B$, defines
the reduced density matrix
$\rho=\mathrm{tr}_{B}\left|\Psi\right\rangle\\!\\!\left\langle\Psi\right|$.
The entanglement Hamiltonian $H_{\text{E}}=-\ln\rho$ has eigenvalues
$\\{E\\}$, known as the entanglement spectrum of $\rho$ Haldane08 . The
entanglement spectrum captures the correlations in the ground state. Moreover,
the universal properties of the actual Hamiltonian of the interacting system
are reflected in $H_{\text{E}}$ Haldane08 ; Calabrese ; QiKatsuraLudwig ;
QiWuZhang . For this reason, we diagnose the effect of interactions as well as
identify the optimal free model exclusively through ground state correlations.

We introduce the interaction distance between the interacting, $\rho$, and the
free, $\sigma$, reduced density matrices

$\displaystyle D_{\cal F}(\rho)=\min_{\sigma\in{\cal F}}D(\rho,\sigma),$ (1)

where $D(\rho,\sigma)=\frac{1}{2}\mathrm{tr}\sqrt{(\rho-\sigma)^{2}}$ is the
trace distance and ${\cal F}$ is the manifold of all free-fermion states. Note
that unlike previous works Genoni:2008el ; Marian:2013cj ; Gertis:2016wj ;
Schilling ; Vollhardt ; Nonfreeness ; ClosestSlater , the manifold $\cal F$
contains all Gaussian states in any set of fermionic quasiparticle operators
$\\{c\\}$. The trace distance has a physical interpretation in terms of
distinguishability between $\rho$ and $\sigma$ when measuring observables
Fuchs . Alternative state-distance measures can equally well be employed
NielsenChuang . The quantity $D_{\cal F}$ has a geometric interpretation as
the distance of the density matrix $\rho$ from ${\cal F}$.

In order to compute $\mathcal{D}_{\mathcal{F}}$, we need to perform the
minimisation in Eq. (1). We can show that this problem can be reduced to
varying only the spectrum of $\sigma$. Consider the basis where $\rho$ is
diagonal and its entanglement spectrum $\\{E\\}$ is available. A general
$\sigma\in{\cal F}$ need not be diagonal in the same basis as $\rho$. However,
the trace distance between $\rho$ and $\sigma$ is minimised when $\sigma$
commutes with $\rho$, i.e., it is simultaneously diagonal with $\rho$ and
their eigenvalue spectra are rank ordered Markham . Indeed, if there existed a
$\sigma$ which minimised $D_{\cal F}$ but did not commute with $\rho$, then
that minimum would not be a global one. As a consequence, the eigenstates of
the optimal model $\sigma$ are the same as the eigenstates of $\rho$. So the
minimisation for obtaining $D_{\mathcal{F}}$ now involves a variation with
respect to the eigenvalues of $\sigma$ that can be given in terms of its
entanglement spectrum.

Having to minimise only with respect to the spectrum of $\sigma$ is a
significant simplification of the optimisation problem. A further exponential
simplification is possible if we consider the structure of its entanglement
spectrum. Consider $N$ single-particle entanglement energies $\\{\epsilon\\}$
with allowed occupations $\\{n_{i},i=1,..,N;n_{i}=0,1\\}$ corresponding to the
independent modes $\\{c\\}$ obeying Fermi-Dirac statistics. The full
entanglement spectrum $\\{E^{f}\\}$ of $\sigma$ contains exponentially many
levels, $2^{N}$, as a function of system size, $N$. However, a special
property of free systems is that due to Wick’s theorem their entanglement
spectrum can be built from a set of single-particle entanglement energies
$\epsilon_{i}$ Peschel03  according to

$E^{\text{f}}_{k}(\\{\epsilon\\})=E_{0}+\sum_{i=1}^{N}n_{i}(k)\,\epsilon_{i},$
(2)

where $E_{0}$ is a normalisation constant. The index $k$ runs over the many-
body spectrum, and for each $k$ there are associated occupations
$n_{i}(k)\in\\{0,1\\}$. Hence, the interaction distance, $D_{\cal F}$, can be
cast as a minimisation with respect to the $N$-many single-particle energies

$D_{\cal
F}(\rho)=\min\limits_{\\{\epsilon\\}}\frac{1}{2}\sum\limits_{k}\left|e^{-E_{k}}-e^{-E^{\text{f}}_{k}(\\{\epsilon\\})}\right|.$
(3)

As $N$ increases at most linearly with system size, expression (3) provides
the means to efficiently compute the interaction distance and obtain the
optimal free model of any state of an interacting theory whenever its
entanglement spectrum $\\{E\\}$ is accessible.

According to (3) the interaction distance is zero when the entanglement
spectrum of $\rho$ satisfies the combinatorial structure given in (2) for
certain single-particle energies, $\\{\epsilon\\}$. This generalises our
concept of free models. In general, no $\epsilon$’s exist that satisfy all
these constraints as their number grows exponentially with the system size,
while the number of $\epsilon$’s grows only linearly. Due to the properties of
the trace distance NielsenChuang  we have that the interaction distance takes
values $D_{\cal F}\in[0,1]$. The condition $D_{\cal F}=0$ corresponds to a
system that can be exactly described by the free fermions $\\{c\\}$, while
$D_{\cal F}=1$ is the maximal distance a state can be from a free description.
When $D_{\cal F}$ is approximately zero, then the deviation in expectation
values or correlation functions between $\rho$ and the approximation $\sigma$
is bounded. In particular, the interaction distance is sensitive only to
deviations in the low lying entanglement spectrum, which dominate the
expectation values of physical observables.

At this point we might wonder what significance the bipartition holds in
definition (1) of the interaction distance. From a quantum information point
of view, the partial trace serves as a quantum channel through which we view
the ground state. This channel is Gaussian if it maps a Gaussian state to
another Gaussian state RevModPhys.84.621 . The combinatorial structure of (2)
describes the eigenvalue spectrum of any Gaussian density matrix. In matching
the spectrum we are testing compatibility between canonical operators in which
the original state is a Gaussian state and the quantum channel is Gaussian,
and our ability to describe the state in terms of modes separable to $A$ and
$B$. If there exist such modes then $D_{\mathcal{F}}$ is zero.

A useful insight in the form of the modes $\\{c\\}$ is given by a constructive
derivation of (3) from (1). The canonical algebra of $\\{c\\}$ is invariant
under arbitrary unitary transformations on $\sigma$. Hence these
transformations keep $\sigma$ in ${\cal F}$, and $\mathcal{F}$ contains all of
its unitary orbits. Optimising for the minimum $D_{\mathcal{F}}$ involves the
following steps. First, $\mathcal{F}$ is decomposed into equivalence classes
which are the mentioned unitary orbits. Within each class, the trace distance
is minimised by a certain representative $\sigma$ which commutes with $\rho$
Markham . Then, $D_{\mathcal{F}}$ is obtained by taking the minimum over
representatives of each class. Since within these unitary orbits the trace
distance is minimised when $\sigma$ and $\rho$ are simultaneously diagonal,
the free modes $\\{c\\}$ are given as Schmidt vectors corresponding to the
single-particle entanglement levels $\\{\epsilon\\}$ for which we have
optimised. Note that the unitarily transformed modes are, in general, a non-
linear combination of the original modes, though they still describe a free
model. From the optimal state $\sigma$, we can identify an effective free
physical Hamiltonian in terms of the emergent quasiparticles $\\{c\\}$ and we
refer to this as the optimal free model. It can be found via a non-unique
procedure Fidkowski ; HolevoWerner  which utilises the two-point correlations
of $\\{c\\}$ with respect to $\sigma$.

General optimal free description:– To determine the value of $D_{\cal F}$ we
have chosen in (2) the statistics of the free quasiparticles to be fermionic.
Generalising further, we could allow for an optimal free description with
respect to other quasiparticle statistics than the constituent fermions of the
original interacting system. However, systems that comprise quasiparticles
with different mutual statistics have in general Hilbert spaces of unequal
dimension. This can be directly resolved in the following way. The
entanglement spectrum of a pure state $\left|\psi\right\rangle$ can be found
from the Schmidt decomposition, $\left|\psi\right\rangle\cong
V_{A}\Sigma{}V_{B}^{\dagger}$, where the state has been reinterpreted as a map
between the Hilbert spaces of the $A$ and $B$ parts. The isometries $V_{A}$
and $V_{B}$ map from $A$ to $B$ through an intermediate space
$\mathcal{H}_{\text{E}}$, the entanglement Hilbert space on which the matrix
of singular values $\Sigma$ and the related entanglement Hamiltonian,
$H_{\text{E}}$, are defined. One can increase the dimension of
$\mathcal{H}_{\text{E}}$ by the addition of zero singular values to $\Sigma$
and the addition of linearly independent columns to $V_{A}$ and $V_{B}$ so
that they remain full-rank, without changing the correlations in
$\left|\psi\right\rangle$. If the dimension of $\mathcal{H}_{\text{E}}$ is
made to exceed the Hilbert space dimension of either subsystem then the
isometries serve as projectors removing non-physical states with high
entanglement energy. This fact allows us to compare interacting systems with
optimal free models that may have different Hilbert space dimension.

Efficiency of computing $D_{\cal F}$:– In the following we shall demonstrate
with two examples that the interaction distance, $D_{\cal F}$, is a versatile
quantity that can be evaluated numerically or analytically. Importantly, the
calculation of $D_{\mathcal{F}}$ requires only the knowledge of the ground
state, thus it is more efficient than other possible measures e.g., based on
directly comparing the structure of actual energy spectra.

When $\rho$ is the reduced density matrix of a gapped 1D system, then $D_{\cal
F}$ can be numerically determined efficiently with $L$, the linear size of the
system. Gapped 1D states have area-law entanglement Hastings07  which bounds
the number of significant levels in the entanglement spectrum in the
thermodynamic limit. Below we use the well known density matrix
renormalisation group algorithm (DMRG) White  to efficiently obtain ground
states of finite 1D systems. Moreover, in this case the low lying
$\\{\epsilon\\}$ in (2), and as a consequence the state $\sigma$, will
converge exponentially fast with $L$. Thus, the minimisation procedure in (3)
is also efficient (see Methods), as it need only involve a finite number of
significant levels determined by the correlation length, even in the
thermodynamic limit.

For critical 1D systems, logarithmic corrections to the area law are possible,
which leads to the polynomial complexity in determining $\rho$. For critical
1D states, a multi-scale renormalisation ansatz VidalMERA  can be implemented
in order to obtain the entanglement spectrum. If the system is gapless, the
number of significant entanglement levels will increase but only polynomially
with system size $L$, hence the optimisation procedure for determining the
single-particle energies $\epsilon$ remains efficient.

Finally, in higher-dimensional systems, our method is reliant on the
efficiency of the current methods in the literature for computing the
entanglement spectrum of the ground state. For 2D systems, one can use
iterative methods such as the Lanczos algorithm in order to access only the
ground state in the exact diagonalisation framework. Furthermore, Monte-Carlo
algorithms QMC  and 2D tensor networks Verstraete ; Orus  can be used in a
variety of systems to variationally approximate the ground state. Then the
typical runtime complexity of $D_{\mathcal{F}}$ for an entanglement spectrum
from a disk partition is polynomial in the correlation length and thus
efficient.

Finite size scaling of $D_{\cal F}$:– The quantity $D_{\mathcal{F}}$, through
its dependence on the entanglement spectrum, inherits the information about
both short- and long-wavelength properties of the system. As pointed out by Li
and Haldane Haldane08 , the entanglement spectrum of a gapped phase exhibits a
generic separation into the universal long-wavelength part and a non-universal
short-distance part, the two being separated by the entanglement gap Haldane08
. Assuming that the linear size of the system’s quasiparticles, $\ell$, is
much smaller than the linear size of the partition $A$, the long-wavelength
information corresponds to correlated quasiparticle excitations across the
entanglement partition, while the short-distance physics is associated with
internal structure of the quasiparticles. The non-universal part is thus a
boundary effect which is insensitive to variation in the subsystem size. In
the thermodynamic limit, the non-universal part is exponentially suppressed in
a gapped phase Haldane08 , as seen from (3), and $D_{\mathcal{F}}$ then
predominantly describes the universal properties of the system.

At critical points where the quasiparticles remain well defined, i.e., $\ell$
stays finite, a large but finite system of linear size $L\gg\ell$ chops off
some of the correlations between the quasiparticles. We surmise that the
finite-size scaling of $D_{\mathcal{F}}$ near such critical regions follows
the ansatz

$D_{\cal
F}\approx{(L^{-1}+\theta)}^{\zeta}f\big{(}(g-g_{\text{c}})L^{1/\nu}\big{)},$
(4)

where $f$ is an undetermined function, and $\zeta$ and $\nu$ are the critical
exponents. The constant $\theta\geq 0$, which vanishes in the standard power-
law scaling ansatz FisherBarber , accommodates the fact that the interaction
distance is bounded from above, $D_{\mathcal{F}}\in[0,1]$. A simple scaling
analysis (see Methods) shows that $\nu$ is the correlation length exponent,
while $\zeta$ determines the effect of interactions in the renormalisation
group sense. For example, when interactions are relevant, $D_{\cal F}$ should
remain non-zero as $L$ increases, which dictates $\zeta\leq 0$. On the other
hand, when interactions are irrelevant, it is expected that $D_{\cal F}$
decreases with $L$ near the critical regions, in which case $\zeta>0$. Note,
however, that it is possible for interactions to be irrelevant and still yield
finite $D_{\mathcal{F}}$ in the thermodynamic limit. This is because,
$D_{\mathcal{F}}$ may be sensitive to non-universal (short distance)
properties of the system, which can give a residual non-zero contribution
parametrised by $\theta$.

Application to Ising model:– For concreteness we consider the example of the
1D ferromagnetic (FM) and antiferromagnetic (AFM) Ising model in both
transverse, $h_{z}$, and longitudinal, $h_{x}$, fields (see Ref. dutta2015
for a recent review). By using exact diagonalisation to determine the
entanglement spectrum, we compute $D_{\mathcal{F}}$ across the phase diagram
and examine its scaling around criticality and its convergence in the
thermodynamic limit. Finally, we find optimal free-fermion models for each
point $(h_{z},h_{x})$ of the phase space.

Figure 1: Interaction distance of the quantum Ising model. (a) Interaction
distance across the phase diagram of the ferromagnetic and (b)
antiferromagnetic models at system size $L=16$ with periodic boundary
conditions. The interaction distance takes non-negligible values only in the
vicinity of the critical point and critical line which are sketched in red.
The scaling behaviour of $D_{\cal F}$ along the dashed lines is given in Fig.
2.

The interacting Hamiltonians are given by

$H_{\pm}=-\sum_{j=1}^{L}(\pm\sigma_{j}^{x}\sigma_{j+1}^{x}+h_{z}\sigma_{j}^{z}+h_{x}\sigma_{j}^{x}),$
(5)

where $H_{+}$ stands for FM and $H_{-}$ for AFM with periodic boundary
conditions. In the presence of only transverse field ($h_{x}=0$), model (5)
maps to free fermions via the Jordan-Wigner transformation. A non-zero
longitudinal field, $h_{x}$, introduces non-local interactions between
fermions. A quantum critical point at $h_{z}=1$ separates an ordered and a
disordered phase of the free system which are related by a self-duality
KramersWannier . The FM model has a single critical point at
$(h_{z}=1,h_{x}=0)$, while the AFM model has a critical line connecting
$(h_{z}=1,h_{x}=0)$ with the classical point $(h_{z}=0,h_{x}=2)$
Ovchinnikov2003 .

Figure 2: Finite-size scaling of the interaction distance for the Ising model.
$D_{\mathcal{F}}$ for a number of systems sizes of (a) the ferromagnetic model
and (c) the antiferromagnetic model along the dashed lines given in Fig. 1. We
obtain the critical exponents $\nu_{\text{FM}}\approx 0.533$,
$\zeta_{\text{FM}}\approx-1.4$ and $\nu_{\text{AFM}}\approx 1.252$,
$\zeta_{\text{AFM}}\approx 1.11$ by fitting to (4). (b) and (d) show the
scaling collapse of $D_{\cal F}$, where $x_{L}=(g-g_{\text{c}})L^{1/\nu}$ and
$y_{L}=D_{\mathcal{F}}L^{\zeta}$. This demonstrates the applicability of the
scaling ansatz (4).

Minimising the interaction distance over the phase diagram we find that
$D_{\mathcal{F}}$ decays with $L$ away from critical regions as shown in Fig.
1, with the variational parameters $\\{\epsilon\\}$ converging exponentially
(see Methods). Thus the model can be faithfully described by a free theory in
these regions of the phase diagram. The exceptions only occur infinitesimally
close to the FM critical point and at the AFM classical critical point. This
is remarkable because these models are non-integrable and a priori have strong
quantum fluctuations due to all energy scales being comparable in magnitude.

Strong correlations build up near criticality where the effect of interactions
is most significant and $D_{\cal F}$ can take higher values. These regions,
however, shrink around criticality as $L$ increases. To examine this we
consider the finite-size scaling of $D_{\cal F}$ using ansatz (4), as shown in
Fig. 2, along the paths $h_{z}=1$ (FM) and $h_{x}=h_{z}$ (AFM) shown in Fig.
1. Note that, despite being near criticality, the values of
$\mathcal{D}_{\mathcal{F}}$ remain small ($\mathcal{D}_{\mathcal{F}}\ll 1$)
for both models, thus we set $\theta=0$ in the ansatz (4). We obtain critical
exponents $\zeta_{\text{FM}}\approx-1.4$ and $\zeta_{\text{AFM}}\approx 1.11$,
showing that interactions have a dramatically different effect in the two
models. In the FM case the interactions are a singular perturbation to the
critical point, whereas in the AFM case their effect diminishes as $L$
increases. This behaviour is consistent with the shrinking of the significant
high-$D_{\mathcal{F}}$ regions around criticality. Furthermore, the critical
exponent $\nu_{\text{FM}}\approx 0.533$, is approximately equal to that of the
correlation length $\nu^{\xi}_{\text{FM}}=8/15$ Zamolodchikov . Similarly
$\nu_{\text{AFM}}\approx 1.252$ which is within $20\%$ accuracy to the
correlation length critical exponent for the same cut
$\nu^{\xi}_{\text{AFM}}\approx 1.052$, which we obtain numerically. We account
for the deviation in $\nu_{\text{AFM}}$ by small system sizes and the fact
that we are not perturbing with an operator that has a well-defined scaling
dimension. The critical scaling behaviour of $D_{\mathcal{F}}$ is
independently verified by employing the variational DMRG method rather than
exact diagonalisation, that can efficiently probe significantly larger system
sizes (see Methods).

Figure 3: Mapping the optimal free model to a transverse-field Ising model.
Panel (a) is for the ferromagnetic model and panel (b) for the
antiferromagnetic model both at system size $L=16$. Contours indicate the
transverse-field $h_{z}^{\text{f}}$, dashed lines for those in the symmetry
broken phase and solid otherwise (including the critical value
$h_{z}^{\text{f}}=1$). The background colour plot gives the distance
$\min_{h^{\text{f}}_{z}}D(\sigma,\sigma^{\text{f}})$ (log scale) signifying
the success of the mapping. In this way the interacting system is given a
description in terms of a free-fermion model. The region near the $h_{z}=0$
classical axis of the AFM is removed (grey region) because our calculations do
not resolve all symmetries of the system.

Finally, we are in position to identify the optimal free model that describes
the interacting system given by an instance of (5). In particular, we identify
the free Ising model $H_{\pm}(h_{z}^{\text{f}},0)$ with transverse field
$h_{z}^{\text{f}}$, whose ground state’s entanglement spectrum matches
$\sigma$’s obtained from (1) for each point $(h_{z},h_{x})$. This is simply
obtained by minimising
$D\big{(}\sigma(h_{z},h_{x}),\sigma^{\text{f}}(h_{z}^{\text{f}},0)\big{)}$
over $h_{z}^{\text{f}}$. As a result we observe that in the FM case, adding
infinitesimal interactions to the free Ising model with $h_{z}<1$ maps the
model to a free Ising with $h_{z}>1$ in a discontinuous way, as shown in Fig.
3 (a). When $h_{z}>1$, the introduction of interactions maps the model to a
neighbouring free model continuously. In the AFM case the interactions are
irrelevant. Indeed, Fig. 3 (b) shows that the whole phase diagram maps
trivially to the free model even very near criticality.

The distance $\min_{h^{\text{f}}_{z}}D(\sigma,\sigma^{\text{f}})$ is shown by
the colour scale in Fig. 3. We see that away from criticality $\sigma$ can be
mapped to the free Ising model with a high fidelity. In the thermodynamic
limit we expect the critical line of the AFM to also be identified with a free
Ising model because $\min_{h^{\text{f}}_{z}}D(\sigma,\sigma^{\text{f}})$
decreases with $L$ and the conformal field theory, which describes the point
$(h_{z}=1,h_{x}=0)$, also governs the entire critical line dutta2015 . This
analysis reveals that the optimal free model is local and fermionic throughout
the phase diagrams, except at the critical point of the FM model and the
classical critical point of the AFM model whose ground state is
macroscopically degenerate.

Maximally interacting model:– In the above analysis, we found for the Ising
model that the interaction distance vanishes in the thermodynamic limit in
almost the entire phase diagram. For that model we identified the optimal free
fermion model that effectively describes it throughout the phase space. We now
present a truly interacting model that cannot be approximated by free
fermions, i.e., has a non-zero $D_{\mathcal{F}}$ in the thermodynamic limit,
even away from criticality. In our analytical approach, we first construct the
entanglement spectrum that gives a non-zero $D_{\mathcal{F}}$ with respect to
free fermions. Then, we identify the physical system and emergent
quasiparticles that correspond to the derived entanglement spectrum.

We consider the simple case where the entanglement Hamiltonian comprises two
fermionic modes. We aim to maximise the interaction distance with respect to
the entanglement spectrum $\\{E_{\text{max}}\\}$ which consists of four
levels. We perform a maximisation of the interaction distance from free-
fermion entanglement spectra generated by two single-particle energies,
$\epsilon_{1},\epsilon_{2}$. It can be analytically proven (see Supplementary
Note 2) that the maximum interaction distance in this case is

$D_{\mathcal{F}}(\rho_{\text{max}})=\max_{\\{\lambda\\}}{D_{\mathcal{F}}}(\rho)=\frac{1}{6},$
(6)

where the maximisation is performed with respect to the eigenvalues
$\\{\lambda\\}$ of $\rho$. The spectrum of the reduced state
$\rho_{\text{max}}$ that maximises $D_{\mathcal{F}}(\rho)$ is the maximally
mixed rank-$3$ spectrum
$\\{\lambda_{\text{max}}\\}=\\{\frac{1}{3},\frac{1}{3},\frac{1}{3},0\\}$, with
entanglement spectrum $\\{E_{\text{max}}\\}=\\{\ln 3,\ln 3,\ln 3,\infty\\}$.

We would now like to find the parent physical system which exhibits such
correlations in the ground state so that it saturates the maximum of
$\mathcal{D}_{\mathcal{F}}$. Since the entanglement spectrum
$\\{E_{\text{max}}\\}$ has a three-fold degeneracy, it is natural to consider
models which support fractionalised excitations. In particular, a
$\mathbb{Z}_{3}$ quantum clock model effectively describes the edge physics of
a 2D topological phase Fendley  and can in principle be realised in the
laboratory Fendley ; Klinovaja ; Vaezi ; Cheng ; YouWen . The $\mathbb{Z}_{3}$
Hamiltonian in the topological phase at its renormalisation fixed point is

$H_{\mathbb{Z}_{3}}=-\sum_{j}{\tau_{j}^{\dagger}\tau_{j+1}+\mathrm{h.c.}},$
(7)

where $\tau_{j}$ are non-Hermitian clock operators which commute between sites
and can be represented locally as
$\tau_{j}=\mathrm{diag}\left(1,e^{i2\pi/3},e^{-i2\pi/3}\right)$ satisfying
$\tau_{j}^{3}=\mathbb{1}$.

The ground state of (7) hosts topologically protected parafermionic zero-modes
exponentially localised on open boundaries ParaFendley . These correspond to
three-fold degeneracy in the entanglement spectrum
$\\{E\\}=\\{\ln{3},\ln{3},\ln{3}\\}$ from maximally entangled pairs of
parafermions across each virtual boundary. We have seen in a previous section
that we can increase the dimension of the entanglement Hilbert space ${\cal
H}_{\text{E}}$ by introducing completely uncorrelated states that correspond
to infinite entanglement energy. Hence, the entanglement spectrum of the
$\mathbb{Z}_{3}$ model reproduces the aforementioned $\\{E_{\text{max}}\\}$.
As a result we have readily identified a truly interacting model that gives
$D_{\mathcal{F}}(\rho_{\text{max}})={1/6}$ with respect to free fermions in
the thermodynamic limit. Importantly, we have managed to analytically identify
the ground state of the corresponding interacting model by considering only
the structure of its correlations.

##  II Discussion

In this paper we have identified, for a given interacting theory, the optimal
free state which is closest to its ground state and introduced the interaction
distance between them. Our approach extends beyond mean-field theory, which is
valid for weak correlations and spatial dimensions above the upper-critical
dimension, or perturbation theory which requires weak couplings. Optimising
with respect to the correlations present in ground state captures faithfully
the low energy physics of the interacting system reproducing the observables
with bounded error.

We have numerically determined the interaction distance for the 1D Ising model
in the presence of transverse and longitudinal fields. We used both exact
diagonalisation as well as DMRG, thus demonstrating the efficiency in
determining $D_{\cal F}$ for large 1D gapped systems using standard
techniques. Further, our diagnostic $D_{\mathcal{F}}$ shows that the ground
state of this interacting model in its gapped phases is well described by a
free state, requiring exponentially less information to represent than the
exact ground state. We expect that this could be generalised to an efficient
algorithm for finding a representative (nearly) free state, similar to the
DMRG method which constructs a low Schmidt rank approximation to gapped 1D
ground states. Alternatively, our method can be combined with analytical wave
functions, as in the case of the Bethe ansatz Bethe  or the trial wave
functions in the quantum Hall effect Laughlin ; ReadRezayi ; MooreRead .

We have also verified that there exist quantum states for which $D_{\cal
F}\neq 0$, such as the $\mathbb{Z}_{3}$ quantum clock model. Other possible
candidates are systems that give rise to exotic phases like
high-$T_{\text{c}}$ superconductors and states with intrinsic topological
order, where interactions play a crucial role. This could establish $D_{\cal
F}$ as an interaction order parameter identifying truly interacting systems,
in terms of fermions or bosons, from nearly free ones like the Ising model.
Furthermore, as we have shown for the Ising Model and the $\mathbb{Z}_{3}$
model, it is possible to directly identify a parent Hamiltonian for the
optimal free state.

In introducing our interaction distance, we have generalised the meaning of
freedom in many-body quantum states by choosing the mutual statistics of the
free modes we are optimising over. Varying the statistics of the free modes
used in our optimisation corresponds to changing the free manifold ${\cal F}$
from which the interaction distance is measured. Our construction can be
immediately generalised to soft-core bosons by promoting the occupation $m$ of
each mode in Eq. (2) to a variational parameter taking values $1\leq
m<\infty$. Taking this further, the single-body levels themselves can become
occupation-dependent. Such a modification could accommodate fractionalised
excitations in strongly-correlated states Simon . Another interesting
generalisation would be to introduce the notion of entanglement temperature
Chandran  which shifts the sensitivity of the measure to other parts of the
entanglement spectrum, reflecting entanglement on different length scales.
Finally, we mention that methods for measuring entanglement spectra in optical
lattices have recently been proposed Pichler . Hence, the interaction distance
can be determined for exotic states realised in cold atom systems.

##  III Methods

Optimisation:–  The optimisation to find $\sigma$ and $D_{\mathcal{F}}$ in (3)
is performed by a Monte Carlo basin-hopping strategy Wales1997  using the
Nelder-Mead simplex algorithm for local minimisation within basins of the cost
function $D_{\mathcal{F}}$. This global strategy was selected to counteract an
observed tendency for local methods to get trapped in local minima. The
initial guess for this search is found as follows. The normalisation constant
$E_{0}$ is the lowest entanglement energy of the input entanglement spectrum
$\\{E\\}$. We iteratively construct an approximate set of single-particle
entanglement energies starting from an empty set. First we take the lowest
remaining level in the spectrum and subtract $E_{0}$ to produce a new single-
particle level $\epsilon_{k}$. Then we remove the many-body levels which are
closest to the new combinatorial levels generated according to (2) by the
additional single-particle level. This process is repeated until the input
spectrum is exhausted. We can also introduce a truncation of the entanglement
spectrum cutting off high entanglement energies, making the construction of
the initial guess terminate faster. The minimisation of
$D(\sigma,\sigma^{\text{f}})$ in order to identify the optimal free model for
the Ising Hamiltonian (5) is calculated using a local Nelder-Mead method.

Finite size scaling:–  We perform the finite-size scaling according to an
ansatz (4). The parameters of the collapse were estimated using the method of
Ref. Collapse . From the scaling ansatz (4) and for a trial set of scaling
parameters $g_{\text{c}}$, $\nu$ and $\zeta$, the scaled values
$x_{L}=(g-g_{\text{c}})L^{1/\nu}$ and
$y_{L}=D_{\mathcal{F}}(1/L+\theta)^{-\zeta}$ are calculated from each unscaled
data point $(g,D_{\mathcal{F}})$. From this collection of scaled data points
${(x_{L},y_{L})}$ across all $L$ we implicitly define a so-called master curve
which best represents them. This curve $y(x)$ is defined around a point $x$ as
the linear regression calculated by taking the scaled data points immediately
left and right of $x$ for each system size $L$. We characterise the deviation
of the scaled points $(x_{L},y_{L})$ from the master curve $y(x_{L})$ using
the chi-square statistic $\chi^{2}$. This measure is used as the cost function
for an optimisation problem over the scaling parameters $g_{\text{c}}$, $\nu$,
$\zeta$ and $\theta$ which can be solved using the same techniques as the
previous problems.

Figure 4:  Interaction distance $D_{\mathcal{F}}$ calculated from DMRG. (a)
Density matrix renormalisation group results retaining $\chi$ states for
$\chi=32$ (dots), $64$ (lines) and $128$ (crosses) states for the
ferromagnetic Ising model. The result is already converged for $\chi=64$ for
the biggest system size $L=1024$. (b) Scaling collapse according to ansatz (4)
with the critical exponents $\nu\approx 0.533$, $\zeta\approx-1.4$ and
saturation parameter $\theta\approx 0.025$. The collapsed variables
$x_{L}=(g-g_{\text{c}})L^{1/\nu}$ and
$y_{L}=D_{\mathcal{F}}(1/L+\theta)^{-\zeta}$. (c) Comparison between
$D_{\mathcal{F}}$ calculated for entanglement spectra provided via DMRG
(points), for $\chi=16$ entanglement levels, and exact diagonalisation (lines)
for the same system. The results for these system sizes accessible by exact
diagonalisation agree. All results in this figure are calculated for open
boundary conditions where DMRG is most natural.

In Fig. 4 (a) we show $D_{\mathcal{F}}$ calculated from the entanglement
spectrum using DMRG White , which extends our results from Fig. 2 in the
ferromagnetic case to larger system sizes inaccessible to exact
diagonalisation. In Fig. 4 (b) we show the scaling collapse. We obtain
critical exponents $\nu$, $\zeta$ which are consistent with those found with
exact diagonalisation and $\nu$ is consistent with known results of CFT
Zamolodchikov . There are two differences compared to results in Fig. 2.
First, with our larger sizes we become sensitive to the upper bound of
$D_{\mathcal{F}}$, which gives us $\theta\approx 0.025$. That $\theta$ is non-
zero is readily apparent in the saturation behaviour visible in the unscaled
data and is demanded for consistency with an upper bound. Second, in
accordance with common practice, we perform DMRG using open boundary
conditions, which changes the non-universal function $f$ in the scaling ansatz
(4).

To verify our DMRG results we compare $D_{\mathcal{F}}$ obtained by DMRG and
exact diagonalisation with open boundary conditions. We confirm that they are
in excellent agreement, as shown in Fig.4 (c). Between iterations in DMRG, the
ground state is approximated by retaining a reduced number of entanglement
states. They correspond to the greatest Schmidt weights. The interaction
distance is insensitive to these approximations inherent in DMRG because the
induced error is comparable to the truncation error which is typically kept to
be $10^{-14}$. We demonstrate that our DMRG calculations are converged in the
number of retained states $\chi$ in Fig. 4 (a). Close to criticality, due to
the logarithmic growth in entanglement, we need to retain more states for
larger $L$. The result is already converged for $\chi=64$ at the largest $L$
studied, which justifies our approximation in retaining this number of states
in the finite-size scaling. Note that the single-particle entanglement
energies converge exponentially away from criticality exponentially (as
discussed in the Supplementary Material).

##  IV Supplementary Material

In this Supplementary Material we demonstrate the convergence of the optimal
free description with system size and present a proof that the maximum of
interaction distance, among all entanglement spectra of rank 4, is equal to
1/6.

Convergence of Ising single body levels:– Regarding the rate of convergence to
the optimal free description with $L$, we find that the single particle
energies $\\{\epsilon\\}$ converge exponentially fast to their asymptotic
values, as shown in Fig. 5. This is the case even near criticality due to the
finite size induced gap. In the thermodynamic limit, however, the entanglement
spectrum is gapless at criticality Calabrese . Hence the convergence to the
optimal model is expected to be polynomial in $L$. We observe that a power law
convergence can be well fitted around criticality, with the goodness of the
fit decreasing away from it (data not shown).

Figure 5: Convergence of the lowest lying single particle entanglement
energies. Graph (a) shows the single-particle entanglement energies
$\epsilon_{i}$ of the ferromagnetic model and graph (b) for the
antiferromagnetic model as a function of $L$ at the non-critical generic point
$(h_{z}=0.88,\,h_{x}=0.16)$. All low lying energies are exponentially
converging with $L$ to their asymptotic values.

Maximum interaction distance among rank-4 entanglement spectra:– Consider a
general 4-level entanglement spectrum of a normalised $\rho$ which is
expressed in terms of the probabilities
$\\{p\\}=\\{p_{1},p_{2},p_{3},p_{4}\\}$. The probabilities are ordered to
satisfy the constraints

$\displaystyle p_{1}\geq p_{2}\geq p_{3}\geq p_{4},$ $\displaystyle
p_{1}+p_{2}+p_{3}+p_{4}=1.$ (8)

We begin by parametrising the probabilities corresponding to 2-mode
(equivalently rank-4) free entanglement spectra by $z,a,b$ with one parameter
(we here pick $b$) fixed by normalisation.

Because the probability spectrum has a definite ordering $z\geq za\geq zb\geq
zab$ for any $z,a$, the cost function can be written in a simple algebraic
form,

$\displaystyle 2D(\\{z,a\\},\\{p\\})=$
$\displaystyle\left|z-p_{1}\right|+\left|za-p_{2}\right|$
$\displaystyle+\left|zb-p_{3}\right|+\left|zab-p_{4}\right|.$ (9)

The interaction distance is a minimisation of this cost function,

$D_{\mathcal{F}}(\rho)=\min\limits_{z,a}D(\\{z,a\\},\\{p\\}).$ (10)

We can fix an element of the variational class $\mathcal{F}$ by choosing
$z=p_{1}$ and $a=p_{2}/p_{1}$ which then forms an upper bounding surface on
$D_{\mathcal{F}}(\rho)$ over $\\{p\\}$,

$D_{\text{B}}(\rho)=D(\\{p_{1},p_{2}/p_{1}\\},\\{p\\})\geq
D_{\mathcal{F}}(\rho).$ (11)

Using the normalisation constraint

$b=\frac{1}{z(1+a)}-1=\frac{1}{p_{1}+p_{2}}-1,$ (12)

we can simplify the upper bound surface to

$\displaystyle 2D_{\text{B}}(\rho)$ $\displaystyle=\left|zb-
p_{3}\right|+\left|zab-p_{4}\right|$ (13)
$\displaystyle=2\left|\frac{p_{1}}{p_{1}+p_{2}}-p_{1}-p_{3}\right|$ (14)

We are now interested in finding the maximum of $D_{\text{B}}(\rho)$ with
respect to $\\{p\\}$ with a view to bounding the maximum of
$D_{\mathcal{F}}(p)$. We simplify this problem by instead considering the
square $D_{\text{B}}^{2}(\rho)$ which is easier to manipulate. The square is a
monotone increasing function, therefore we can equivalently maximise
$D_{\text{B}}^{2}(\rho)$ to find the maximum of $D_{\text{B}}(\rho)$.

To solve this constrained maximisation, we take derivatives first with respect
to $p_{3}$,

$\left(\frac{\partial D_{\text{B}}^{2}}{\partial
p_{3}}\right)_{p_{1},p_{2}}=-2\left(\frac{p_{1}}{p_{1}+p_{2}}-p_{1}-p_{3}\right).$
(15)

The only extremal point of the unconstrained problem is a minimum, therefore
when $D_{B}$ is maximised $p_{3}$ must saturate its constraints. Since the
lower bound $p_{3}=0$ is that of a rank-$2$ spectrum which has
$D_{\mathcal{F}}=0$ the maximum is found for $p_{3}=p_{2}$,

$\displaystyle\left(\frac{\partial D_{\text{B}}^{2}}{\partial
p_{2}}\right)_{\begin{subarray}{c}p_{2}=p_{3},\\\ p_{3}\end{subarray}}$
$\displaystyle=$ $\displaystyle-2\left(1+\frac{1}{(p_{1}+p_{2})^{2}}\right)$
(16) $\displaystyle\times\left(\frac{p_{1}}{p_{1}+p_{2}}-p_{1}-p_{2}\right).$

This derivative vanishes for $p_{1}/(p_{1}+p_{2})-p_{1}-p_{2}$ for which
$D_{\text{B}}=0$ and therefore minimal. Hence to maximise $D_{\text{B}}$,
$p_{2}$ must saturate its constraints. At its lower bound $p_{2}=0$ we find
$D_{\text{B}}(p)=0$ which cannot be the maximum, hence the maximum is found
for $p_{2}=p_{1}$,

$\left(\frac{\partial D_{\text{B}}^{2}}{\partial
p_{1}}\right)_{p_{1}=p_{2}=p_{3}}=-4\left(\frac{1}{2}-2p_{1}\right)^{2}.$ (17)

Because once again the only stationary point is a minimum, the constraints on
$p_{1}$ must be saturated. The lower bound on $p_{1}$ is $1/4$ and the upper
bound is $1/3$. Because $p_{1}=1/4$ gives $D_{\mathcal{F}}(p)=0$, the maximum
is found for $p_{1}=p_{2}=p_{3}=1/3$ and by normalisation $p_{4}=0$.

Substituting these values for $\\{p\\}$ in (11) we have
$D_{\text{B}}(1/3,1/3,1/3,0)=1/6$. Then by direct analytic calculation of
$D_{\mathcal{F}}(\\{1/3,1/3,1/3,0\\})$ one can show that this upper bound is
attained. The conclusion is that amongst rank-$4$ probability spectra
$\\{1/3,1/3,1/3,0\\}$ is the unique maximum of $D_{\mathcal{F}}$ achieving
$D_{\mathcal{F}}=1/6$.

##  V Data Availability

Data sharing is not applicable to this article as no datasets were generated
or analysed during the current study.

##  VI Acknowledgements

We would like to thank M. Barkeshli, B. Doyon, P. Fendley, G. Giedke, J.
Garcia-Ripoll, S. Iblisdir, A. Läuchli, T. Neupert, D. Poilblanc, A.
Polkovnikov, K. Shtengel and S. Simon, for inspiring conversations and useful
comments. This work was supported by the EPSRC grants EP/I038683/1,
EP/M50807X/1 and EP/P009409/1. Statement of compliance with EPSRC policy
framework on research data: This publication is theoretical work that does not
require supporting research data.

##  VII Author Contributions

All authors contributed to developing the ideas, analysing the results and
writing the manuscript. C.J.T. implemented the algorithm.

##  VIII Competing financial interests

The authors declare no competing financial interests.

##  IX Corresponding authors

Correspondence to: Christopher J. Turner and Konstantinos Meichanetzidis.

## References

  * (1) Sutherland, B.  _Beautiful Models: 70 Years of Exactly Solved Quantum Many-body Problems_ (World Scientific, 2004).
  * (2) Feynman, R. P.  Atomic theory of the two-fluid model of liquid helium.  _Phys. Rev._ 94, 262–277 (1954).
  * (3) Bardeen, J., Cooper, L. N. & Schrieffer, J. R.  Theory of superconductivity.  _Phys. Rev._ 108, 1175–1204 (1957).
  * (4) Gutzwiller, M. C.  Effect of correlation on the ferromagnetism of transition metals.  _Phys. Rev. Lett._ 10, 159–162 (1963).
  * (5) Laughlin, R. B.  Anomalous quantum Hall effect: An incompressible quantum fluid with fractionally charged excitations.  _Phys. Rev. Lett._ 50, 1395–1398 (1983).
  * (6) Li, H. & Haldane, F. D. M.  Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum Hall effect states.  _Phys. Rev. Lett._ 101, 010504 (2008).
  * (7) Calabrese, P. & Lefevre, A.  Entanglement spectrum in one-dimensional systems.  _Phys. Rev. A_ 78, 032329 (2008).
  * (8) Qi, X.-L., Katsura, H. & Ludwig, A. W. W.  General relationship between the entanglement spectrum and the edge state spectrum of topological quantum states.  _Phys. Rev. Lett._ 108, 196402 (2012).
  * (9) Qi, X.-L., Wu, Y.-S. & Zhang, S.-C.  General theorem relating the bulk topological number to edge states in two-dimensional insulators.  _Phys. Rev. B_ 74, 045125 (2006).
  * (10) Genoni, M. G., Paris, M. G. A. & Banaszek, K.  Quantifying the non-Gaussian character of a quantum state by quantum relative entropy.  _Phys. Rev. A_ 78, 060303 (2008).
  * (11) Marian, P. & Marian, T. A.  Relative entropy is an exact measure of non-Gaussianity.  _Phys. Rev. A_ 88, 012322 (2013).
  * (12) Gertis, J., Friesdorf, M., Riofrío, C. A. & Eisert, J.  Estimating strong correlations in optical lattices.  _Phys. Rev. A_ 94, 053628 (2016).
  * (13) Schilling, C., Gross, D. & Christandl, M.  Pinning of fermionic occupation numbers.  _Phys. Rev. Lett._ 110, 040404 (2013).
  * (14) Byczuk, K., Kuneš, J., Hofstetter, W. & Vollhardt, D.  Quantification of correlations in quantum many-particle systems.  _Phys. Rev. Lett._ 108, 087004 (2012).
  * (15) Gottlieb, A. D. & Mauser, N. J.  Nonfreeness and related functionals for measuring correlation in many-fermion states (2015).  Preprint at https://arxiv.org/abs/1510.04573.
  * (16) Zhang, J. M. & Kollar, M.  Optimal multiconfiguration approximation of an $n$-fermion wave function.  _Phys. Rev. A_ 89, 012504 (2014).
  * (17) Fuchs, C. A. & van de Graaf, J.  Cryptographic distinguishability measures for quantum-mechanical states.  _IEEE Trans. Inf. Theory_ 45, 1216–1227 (1999).
  * (18) Nielsen, M. A. & Chuang, I. L.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, New York, NY, USA, 2011), 10th edn.
  * (19) Markham, D., Miszczak, J. A., Puchała, Z. & Życzkowski, K.  Quantum state discrimination: A geometric approach.  _Phys. Rev. A_ 77, 042111 (2008).
  * (20) Peschel, I.  Calculation of reduced density matrices from correlation functions.  _J. Phys. A_ 36, L205 (2003).
  * (21) Weedbrook, C. _et al._ Gaussian quantum information.  _Rev. Mod. Phys._ 84, 621–669 (2012).
  * (22) Fidkowski, L.  Entanglement spectrum of topological insulators and superconductors.  _Phys. Rev. Lett._ 104, 130502 (2010).
  * (23) Holevo, A. S. & Werner, R. F.  Evaluating capacities of bosonic gaussian channels.  _Phys. Rev. A_ 63, 032312 (2001).
  * (24) Hastings, M. B.  An area law for one-dimensional quantum systems.  _J. Stat. Mech. Theor. Exp._ 2007, P08024 (2007).
  * (25) White, S. R.  Density matrix formulation for quantum renormalization groups.  _Phys. Rev. Lett._ 69, 2863–2866 (1992).
  * (26) Vidal, G.  Class of quantum many-body states that can be efficiently simulated.  _Phys. Rev. Lett._ 101, 110501 (2008).
  * (27) Zhang, S.  _Quantum Monte Carlo Methods for Strongly Correlated Electron Systems_ (Springer New York, New York, NY, 2004).
  * (28) Verstraete, F., Murg, V. & Cirac, J. I.  Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems.  _Advances in Physics_ 57, 143–224 (2008).
  * (29) Orús, R.  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  _Annals of Physics_ 349, 117 – 158 (2014).
  * (30) Fisher, M. E. & Barber, M. N.  Scaling theory for finite-size effects in the critical region.  _Phys. Rev. Lett._ 28, 1516–1519 (1972).
  * (31) Dutta, A. _et al._ _Quantum Phase Transitions in Transverse Field Models_ (Cambridge University Press, 2015).
  * (32) Kramers, H. A. & Wannier, G. H.  Statistics of the two-dimensional ferromagnet. Part I.  _Phys. Rev._ 60, 252–262 (1941).
  * (33) Ovchinnikov, A. A., Dmitriev, D. V., Krivnov, V. Y. & Cheranovskii, V. O.  Antiferromagnetic Ising chain in a mixed transverse and longitudinal magnetic field.  _Phys. Rev. B_ 68, 214406 (2003).
  * (34) Zamolodchikov, A. B.  Integrals of motion and S-matrix of the (scaled) $T=T_{c}$ Ising model with magnetic field.  _Int. J. Mod. Phys. A_ 04, 4235–4248 (1989).
  * (35) Alicea, J. & Fendley, P.  Topological phases with parafermions: Theory and blueprints.  _Annual Review of Condensed Matter Physics_ 7, 119–139 (2016).  arXiv:eprint http://dx.doi.org/10.1146/annurev-conmatphys-031115-011336.
  * (36) Klinovaja, J., Yacoby, A. & Loss, D.  Kramers pairs of majorana fermions and parafermions in fractional topological insulators.  _Phys. Rev. B_ 90, 155447 (2014).
  * (37) Vaezi, A.  Fractional topological superconductor with fractionalized majorana fermions.  _Phys. Rev. B_ 87, 035132 (2013).
  * (38) Cheng, M.  Superconducting proximity effect on the edge of fractional topological insulators.  _Phys. Rev. B_ 86, 195126 (2012).
  * (39) You, Y.-Z. & Wen, X.-G.  Projective non-abelian statistics of dislocation defects in a ${\mathbb{Z}}_{N}$ rotor model.  _Phys. Rev. B_ 86, 161107 (2012).
  * (40) Jermyn, A. S., Mong, R. S. K., Alicea, J. & Fendley, P.  Stability of zero modes in parafermion chains.  _Phys. Rev. B_ 90, 165106 (2014).
  * (41) Bethe, H.  Zur Theorie der Metalle.  _Zeitschrift für Physik_ 71, 205–226 (1931).
  * (42) Read, N. & Rezayi, E.  Beyond paired quantum Hall states: Parafermions and incompressible states in the first excited Landau level.  _Phys. Rev. B_ 59, 8084–8092 (1999).
  * (43) Moore, G. & Read, N.  Nonabelions in the fractional quantum Hall effect.  _Nuclear Physics B_ 360, 362 – 396 (1991).
  * (44) Davenport, S. C., Rodríguez, I. D., Slingerland, J. K. & Simon, S. H.  Composite fermion model for entanglement spectrum of fractional quantum Hall states.  _Phys. Rev. B_ 92, 115155 (2015).
  * (45) Chandran, A., Khemani, V. & Sondhi, S. L.  How universal is the entanglement spectrum?  _Phys. Rev. Lett._ 113, 060501 (2014).
  * (46) Pichler, H., Zhu, G., Seif, A., Zoller, P. & Hafezi, M.  Measurement protocol for the entanglement spectrum of cold atoms.  _Phys. Rev. X_ 6, 041033 (2016).
  * (47) Wales, D. J. & Doye, J. P. K.  Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms.  _J. Phys. Chem. A_ 101, 5111–5116 (1997).
  * (48) Houdayer, J. & Hartmann, A. K.  Low-temperature behavior of two-dimensional gaussian Ising spin glasses.  _Phys. Rev. B_ 70, 014418 (2004).


supplement.pdf

# Quantum Nonlinear Acoustic Hall Effect and Inverse Acoustic Faraday Effect
in Dirac Insulators

Ying Su<EMAIL_ADDRESS>Center for Integrated Nanotechnology, Los Alamos
National Laboratory, Los Alamos, New Mexico 87545, USA  Alexander V. Balatsky
Nordita, KTH Royal Institute of Technology and Stockholm University 106 91
Stockholm, Sweden  Department of Physics, University of Connecticut, Storrs,
Connecticut 06269, USA  Shi-Zeng Lin<EMAIL_ADDRESS>Theoretical Division, T-4
and CNLS, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA
Center for Integrated Nanotechnology, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA

(August 31, 2024)

###### Abstract

We propose to realize the quantum nonlinear Hall effect and the inverse
Faraday effect through the acoustic wave in a time-reversal invariant but
inversion broken Dirac insulator. We focus on the acoustic frequency much
lower than the Dirac gap such that the interband transition is suppressed and
these effects arise solely from the intrinsic valley-contrasting band
topology. The corresponding acoustoelectric conductivity and magnetoacoustic
susceptibility are both proportional to the quantized valley Chern number and
independent of the quasiparticle lifetime. The linear and nonlinear components
of the longitudinal and transverse topological currents can be tuned by
adjusting the polarization and propagation directions of the surface acoustic
wave. The static magnetization generated by a circularly polarized acoustic
wave scales linearly with the acoustic frequency as well as the strain-induced
charge density. Our results unveil a quantized nonlinear topological
acoustoelectric response of gapped Dirac materials, like hBN and transition-
metal dichalcogenide, paving the way toward room-temperature acoustoelectric
devices due to their large band gaps.

_Introduction._ —Dirac materials, owing to their relativistic energy
dispersion and nontrivial topological properties, have drawn significant
attention [1, 2, 3]. One important feature in many Dirac materials is the
emergence of valley degree of freedom around the low-energy band extrema where
electronic states are well separated in momentum space and form Kramers pairs
in the presence of time-reversal symmetry (TRS) [4, 5]. Similar to spin, the
states with opposite valley indices can be used as binary information carriers
that evokes a surge of interest in valleytronics [6, 7]. Recently, the
acoustic wave has been demonstrated as an efficient method to manipulate the
valley degree of freedom and has led to the valley acoustoelectric (AE) effect
[8, 9, 10], acoustogalvanic effect [11, 12, 13], axial magnetoelectric effect
[14], etc. Within these effects, direct (valley) current and static
magnetization are generated through the nonlinear response to the alternating
strain exerted by different acoustic waves in doped Dirac semimetals and
insulators [8, 9, 10, 11, 12, 13, 14].

The TRS-invariant strain field emerges as a deformation potential as well as a
valley-contrasting gauge field for Dirac fermions [15, 16, 17, 18], that
results in novel quantum effects like the zero-field Landau quantization [18,
19], quantum spin Hall effect [20], valley polarization and inversion [21,
22], etc. Despite the conventional AE effect stemming from the deformation
potential had been identified in 1950s [23, 24], the role of the emergent
gauge field exerted by an acoustic wave passing through the doped Dirac
semimetal or insulator has just been recoginized [11, 12, 13, 14]. In
particular, the interplay between the Dirac fermion and acoustic wave gives
rise to a direct Hall current transverse to the wave vector without the TRS
breaking, similar to the nonlinear Hall effect [25]. To date, the AE Hall
effect arises from the second-order response to the emergent gauge field of
acoustic waves and exists only in the metallic phase [11, 12, 13]. Likewise,
the inverse acoustic faraday effect (IAFE, which generates static
magnetization by a circularly polarized acoustic wave [26, 27]) can be
achieved in Dirac semimetals via the second-order response that results in
rectified magnetization [14].

Discovering the quantum version of various Hall effects is of significant
importance in condensed matter physics by unveiling different topological
phases of insulators and quantized topological responses [28, 29, 30, 31, 32,
33, 34]. However, it remains unclear how the undoped Dirac insulator (i.e.,
the Fermi energy lies in the band gap) responds to different acoustic waves.
In this work, we propose to realize the quantum nonlinear acoustic Hall effect
(QNAHE) and IAFE through the quantized topological AE response in 2D Dirac
insulators with TRS but no inversion symmetry. In QNAHE, we show that a
surface acoustic wave (SAW) passing through the Dirac insulator can generate
both longitudinal and transverse currents. Moreover, their linear and
nonlinear components (up to the second order) are tunable by adjusting the
polarization and propagation directions of the SAW. When the acoustic wave is
circularly polarized and propagating in the out-of-plane direction, there is a
static magnetization in the Dirac insulator subject to inhomogeneous strain,
i.e., the realization of IAFE. Here we focus on the acoustic frequency much
lower than the Dirac gap such that the interband transition is suppressed.
Since no free charge carrier is involved in the AE response, the QNAHE and
IAFE arise solely from the intrinsic valley-contrasting band topology of the
Dirac insulator. In both cases, the corresponding AE conductivity and
magnetoacoustic (MA) susceptibility are proportional to the quantized valley
Chern number and independent of the quasiparticle lifetime, manifesting their
topological origins. Our results provide a new mechanism of rectification that
converts the alternating acoustic wave into direct (Hall) current and static
magnetization in the Dirac insulator via the quantized nonlinear topological
AE response.

Figure 1: (a) Schematic setup for the QNAHE. The upper panel shows the
directions of the SAW vector $\bm{q}$, longitudinal current $j_{\parallel}$,
and transverse current $j_{\perp}$ with respect to the $x$ and $y$ axes that
are fixed along the zigzag and armchair directions of the honeycomb lattice,
respectively. $C_{2\parallel}$ represents the two-fold rotation along the axis
of the SAW wavevector $\bm{q}$. The lower panel depicts a 2D Dirac material
placed on a piezoelectric substrate in which a SAW propagates along its wave
vector direction. (b) and (c) The polar plots of the absolute value of the
linear AE conductivities. (d) The polar plot of the absolute value of the
nonlinear AE conductivities.

_Model and topological acoustoelectric response._ —Here we consider the low-
energy states of a 2D Dirac insulator, like hBN or transition-metal
dichalcogenide, which can be generally described by the Dirac action

$S=\int
d^{3}x\bar{\psi}\left(i\cancel{\partial}+\cancel{A}+\tau_{z}\cancel{a}-\tau_{z}m\right)\psi,$
(1)

where $\cancel{\partial}=\gamma^{\mu}\partial_{\mu}$ denotes the Feynman slash
notation and
$\gamma^{\mu=0,1,2}=\left(\sigma_{z},i\sigma_{y},-i\sigma_{x}\right)$ with
$\sigma_{i}$ being the Pauli matrix for internal degrees of freedom, e.g.,
sublattice or orbital. $\bar{\psi}=\psi^{\dagger}\gamma^{0}$ represents the
Dirac adjoint. For convenience, we set the units $e=\hslash=v_{F}=1$ where
$v_{F}$ is the Dirac velocity. $A_{\mu}=\left(\Phi,-\bm{A}\right)$ is the real
gauge field, while $a_{\mu}=\left(0,-\bm{a}\right)$ is the emergent gauge
field. In the presence of a strain field $\bm{u}(\bm{r})$, the deformation
potential is $\Phi=g_{D}(u_{xx}+u_{yy})$ and the emergent gauge field reads
$\bm{a}=g\left(u_{yy}-u_{xx},2u_{xy}\right)$ where
$u_{ij}=\frac{1}{2}\left(\partial_{i}u_{j}+\partial_{j}u_{i}+\partial_{i}d\partial_{j}d\right)$
is the strain tensor [15, 16, 17]. Here $d(\bm{r})$ is the out-of-plane
displacement field and $g=\beta/2a_{0}$ with the Grüneisen parameter $\beta$
and lattice bond length $a_{0}$. Due to the TRS, the emergent gauge field and
Dirac mass are valley-contrasting, as described by the last two terms of Eq.
(1) where $\tau_{i}$ is the Pauli matrix for valley degrees of freedom.

To investigate how the Dirac insulator responds to the emergent gauge field,
we integrate out the fermion in Eq. (1) which yields a Chern-Simons (CS)
action [35]

$S_{\mathrm{CS}}=\frac{\mathrm{sgn}(m)}{2\pi}\int
d^{3}x\epsilon^{\mu\nu\lambda}A_{\mu}\partial_{\nu}a_{\lambda},$ (2)

where $\mathrm{sgn}(m)=C_{K}-C_{K^{\prime}}$ is the valley Chern number. In
contrast to the conventional CS term that vanishes due to TRS, Eq. (2)
represents a cross CS term between the real and emergent gauge fields, leading
to the topological current density

$j^{\mu}=\frac{\delta S_{\mathrm{CS}}}{\delta
A_{\mu}}=\frac{\text{sgn}(m)}{\pi}\epsilon^{\mu\nu\lambda}\partial_{\nu}a_{\lambda}.$
(3)

In terms of the emergent gauge field, the pseudo-electric and magnetic fields
are defined as $\bm{E}_{s}=-\partial_{t}\bm{a}$ and
$B_{s}=\left(\nabla\times\bm{a}\right)\cdot\hat{\bm{z}}$. Then Eq. (3) can be
re-expressed as

$\bm{j}=\frac{2e^{2}\mathrm{sgn}(m)}{h}\left(\begin{array}[]{cc}0&-1\\\
1&0\end{array}\right)\bm{E}_{s},\quad\rho=-\frac{2e^{2}\mathrm{sgn}(m)}{h}B_{s},$
(4)

where the physical units are recovered. It demonstrates that: (i) there is a
current transverse to the pseudo-electric field; (ii) The charge density is
proportional to the pseudo-magnetic field; (iii) Both current and charge
densities are proportional to the quantized valley Chern number.

In principal, Eqs. (3) and (4) reflect the linear response to the emergent
gauge field and hence the alternating in-plane strain generates only the
alternating linear current in graphene [35, 36]. Meanwhile, the strain-induce
charge density yields the piezoelectricity [37, 38] and magnetoelectric effect
[39] in the Dirac system. Nevertheless, there exists a second-order response
to the strain in Eq. (3) since the emergent gauge field is proportional to the
strain tensor that is quadratic in the gradients of the out-of-plane
displacement $d(\bm{r})$. To the best of our knowledge, this nonlinear
topological AE response stemming from the CS term and acoustic wave in the
Dirac insulator has not been discussed before and is the main focus of this
study by giving rise to the QNAHE and IAFE. Using the realistic material
parameters, we will show that the acoustic-wave-induced direct topological
current and static magnetization in Dirac insulators are measurable with the
magnitudes comparable to other nonlinear responses.

_Quantum nonlinear acoustic Hall effect._ —To trigger the topological current
in Eq. (4), we consider a generic SAW

$\bm{u}=\mathrm{Re}\left[\left(u_{\parallel}\hat{\bm{q}}+u_{\perp}\partial_{\theta}\hat{\bm{q}}-iu_{z}\hat{\bm{z}}\right)e^{i(\bm{q}\cdot\bm{r}-\omega
t)}\right],$ (5)

passing through the Dirac insulator, as shown in Fig. 1(a). Here the 2D
material is placed on top of a piezoelectric substrate and the SAW can be
launched by the interdigital transducer [40]. The wave vector
$\bm{q}=q\left(\cos\theta,\sin\theta\right)$ is measured with respect to the
$x$ and _$y$_ axes that are fixed along the zigzag and armchair directions of
the honeycomb lattice respectively, see Fig. 1(a). When $u_{\parallel}=0$
($u_{\perp}=0$), Eq. (5) describes the circularly polarized (Rayleigh) SAW.
When $u_{\parallel}=u_{z}=0,$ it becomes the Bleustein-Gulyaev SAW [41]. The
topological current driven by the SAW can be decomposed into the longitudinal
and transverse components as
$\bm{j}=j_{\parallel}\hat{\bm{q}}+j_{\perp}\partial_{\theta}\hat{\bm{q}}$ with
$j_{\alpha}=\sum_{n}\mathrm{Re}[j_{\alpha}^{(n)}e^{i(n\bm{q}\cdot\bm{r}-n\omega
t)}]$ whose nonzero components are

$\displaystyle\left(\begin{array}[]{c}j_{\parallel}^{(1)}\\\
j_{\perp}^{(1)}\end{array}\right)=$
$\displaystyle\frac{2e^{2}\mathrm{sgn}(m)gq\omega}{h}\left(\begin{array}[]{c}u_{\parallel}\sin
3\theta+u_{\perp}\cos 3\theta\\\ u_{\parallel}\cos 3\theta-u_{\perp}\sin
3\theta\end{array}\right),$ (10)
$\left(\begin{array}[]{c}j_{\parallel}^{(2)}\\\
j_{\perp}^{(2)}\end{array}\right)=-i\frac{e^{2}\mathrm{sgn}(m)gq^{2}\omega
u_{z}^{2}}{h}\left(\begin{array}[]{c}\sin 3\theta\\\ \cos
3\theta\end{array}\right).$ (11)

A nonlinear (second-order) topological current can be generated by the out-of-
plane component $u_{z}$ of the SAW, besides a linear (first-order) topological
current from the in-plane components $u_{\parallel}$ and $u_{\perp}$.

Now we introduce the AE conductivity tensors for linear and nonlinear currents
as

$j_{\alpha}^{(1)}=\sigma^{\alpha\beta}u_{\beta},\qquad\qquad
j_{\alpha}^{(2)}=\sigma^{\alpha\beta\gamma}u_{\beta}u_{\gamma}.$ (12)

According to Eqs. (10) and (11), the tensors are proportional to the acoustic
frequency as well as the quantized valley Chern number and are independent of
the quasiparticle lifetime. Furthermore, their polar angle dependence exhibits
the $C_{3z}$ symmetry, consistent with the lattice rotational symmetry. In
Figs. 1(b) and 1(c), we show the polar plots of the absolute rank-2 AE
conductivity tensors $|\sigma^{\parallel\alpha}|$ and $|\sigma^{\perp\alpha}|$
(for $\alpha=\parallel$ or $\perp$) that are responsible for the linear
longitudinal current $j_{\parallel}^{(1)}$ and transverse current
$j_{\perp}^{(1)}$ driven by $u_{\alpha}$. Note that the linear currents driven
by $u_{\parallel}$ and $u_{\perp}$ are always perpendicular to each other, as
elucidated in Eq. (10). For the absolute rank-3 conductivity tensors
$|\sigma^{\parallel zz}|$ and $|\sigma^{\perp zz}|$ in Fig. 1(d), they display
the same polar angle dependence as $|\sigma^{\parallel\parallel}|$ and
$|\sigma^{\perp\parallel}|$ in Fig. 1(b). Thus the nonlinear current driven by
$u_{z}$ is always parallel (perpendicular) to the linear current driven by
$u_{\parallel}$ ($u_{\perp}$). In particular for the SAW propagating in the
zigzag (armchair) direction with $\theta=n\pi/3$ ($n\pi/3+\pi/6$),
$\sigma^{\parallel\parallel}=\sigma^{\perp\perp}=\sigma^{\parallel zz}=0$
($\sigma^{\perp\parallel}=\sigma^{\parallel\perp}=\sigma^{\perp zz}=0$). This
is constrained by the symmetry of the system. Under the two-fold rotation
$C_{2\parallel}$ along the axis of the SAW vector $\bm{q}$, as shown in Fig.
1(a), $u_{\parallel}$ and $j_{\parallel}$ ($u_{\perp}$, $u_{z}$ and
$j_{\perp}$) are even (odd). When the SAW propagates in the zigzag (armchair)
direction, the two sublattices are exchanged (invariant) under the
$C_{2\parallel}$ such that the Dirac mass $m$ is odd (even). Because the two
sides of Eqs. (10) and (11) have to transform in the same manner under the
symmetry operations like $C_{2\parallel}$ and TRS, we can explain the
behaviors of the AE conductivities above, see the Supplemental Material (SM)
[42] for detailed symmetry analysis.

Exploiting the anisotropy of AE conductivity tensors, the linear and nonlinear
components of longitudinal and transverse topological currents can be feasibly
tuned by adjusting the polarization and propagation directions of the SAW. For
instance, the Rayleigh (circularly polarized) SAW generates parallel
(perpendicular) linear and nonlinear currents. When the circularly polarized
SAW propagates in the zigzag (armchair) direction, there are linear
(nonlinear) longitudinal current and nonlinear (linear) transverse current. In
practice, as long as the SAW has an out-of-plane component and does not
propagate in the armchair direction, $\sigma^{\perp zz}\neq 0$ leads to a
nonlinear transverse current, i.e., the QNAHE.

The quantized nonlinear topological AE response originating from the quadratic
term in the strain tensor enables not only the frequency multiplier but also
the rectification. To realize the rectified direct topological current, we
consider a pulsed SAW whose amplitude modulation can be encoded in the time-
dependent $u_{\alpha}(t)$. Then the zeroth-order topological current reads

$\left(\begin{array}[]{c}j_{\parallel}^{(0)}\\\
j_{\perp}^{(0)}\end{array}\right)=\frac{e^{2}\mathrm{sgn}(m)gq^{2}u_{z}\partial_{t}u_{z}}{h}\left(\begin{array}[]{c}\sin
3\theta\\\ \cos 3\theta\end{array}\right).$ (13)

By demanding $\partial_{t}\left(u_{z}\partial_{t}u_{z}\right)=0$, we find that
the sawtooth-type SAW with $u_{z}(t)=\bar{u}_{z}\sqrt{(t\mathrm{\;mod}\;T)/T}$
(where $T$ is the period of the pulse) can generate a direct topological
current

$\bm{j}_{\mathrm{DC}}^{(0)}=\frac{e^{2}\mathrm{sgn}(m)gq^{2}\bar{u}_{z}^{2}}{2hT}(\sin
3\theta,\cos 3\theta),$ (14)

within each period and propagating in the same direction as the nonlinear
current in Eq. (11). To estimate the direct current density, we take monolayer
hBN with the Grüneisen parameter $\beta=3.3$ and bond length $a_{0}=1.44$ Å
[37] as an example. For a typical pulsed SAW with $\bar{u}_{z}\sim 10$ nm,
$q\sim 10$ $\mu\mathrm{m}^{-1}$, and $T\sim 1$ $\mu\mathrm{s}$, the current
density is $j_{\mathrm{DC}}^{(0)}\sim 15$ nA/cm, comparable to that in the
doped Dirac insulator [8] and semimetal [11].

Besides the topological current, the emergent gauge field generates a charge
density that satisfies the continuity equation $\partial_{\mu}j^{\mu}=0$, as
demonstrated in Eq. (3). For the SAW in Eq. (5), the charge density is
$\rho=\sum_{n}\mathrm{Re}[\rho^{(n)}e^{i(n\bm{q}\cdot\bm{r}-n\omega t)}]$ with
nonzero components
$\rho^{(1)}=2e^{2}h^{-1}\mathrm{sgn}(m)gq^{2}(u_{\parallel}\sin
3\theta+u_{\perp}\cos 3\theta)$ and
$\rho^{(2)}=-ie^{2}h^{-1}\mathrm{sgn}(m)gu_{z}^{2}q^{3}\sin 3\theta$ according
to Eq. (4). To verify the topological AE response, we also calculate the
charge density from the tight-binding model
$H=\sum_{i}\varepsilon_{i}c_{i}^{\dagger}c_{i}+\sum_{\left\langle
ij\right\rangle}t_{ij}c_{i}^{\dagger}c_{j}$ where $\varepsilon_{i}=\pm m$ is
the sublattice potential and
$t_{ij}=t_{0}\exp\left[-\beta\left(d_{ij}/a_{0}-1\right)\right]$ is the
hopping strength in the honeycomb lattice. Here the bond length $d_{ij}$ is
altered by the SAW and the low-energy states of the lattice model is described
by Eq. (1). Using the realistic parameters of the monolayer hBN, the charge
densities from the CS theory and tight-binding model exhibit excellent
agreement with each other, see SM [42] for details.

_Inverse acoustic Faraday effect._ —The charge density induced by
inhomogeneous strain through the topological AE response generates not only
piezoelectricity [37, 38, 39] but also MA susceptibility. We now evaluate the
static magnetization excited by a circularly polarized acoustic wave in a
strained Dirac insulator, i.e., the realization of IAFE. As displayed in Fig.
2(a), we consider a rippled 2D material on top of a non-flat piezoelectric
substrate and the inhomogeneous strain can be determined by minimizing the
elastic energy functional

$V_{\mathrm{E}}=\int
d^{2}\bm{r}\left[\frac{\lambda}{2}\left(\sum_{i}u_{ii}\right)^{2}+\mu\sum_{i,j}u_{ij}^{2}\right]+\frac{\kappa}{2}\int
d^{2}\bm{r}\left(\nabla^{2}d\right)^{2},$ (15)

where $\lambda$ and $\mu$ are the Lamé factors and $\kappa$ is the bending
rigidity [15, 43]. The topography of the 2D material is encoded in
$d(\bm{r})$. Then the electric charge due to strain-induced emergent gauge
field can be driven into orbital motion by a circularly polarized acoustic
wave
$\bm{u}^{\prime}=\mathrm{Re}[u^{\prime}(\hat{\bm{x}}-i\hat{\bm{y}})e^{i(q_{z}z-\omega
t)}]$ propagating in the out-of-plane direction. This yields a static
magnetization

$\bm{M}_{\bm{k}}=\chi_{\bm{k}}\bm{u}_{q_{z},\omega}^{\prime}\times\bm{u}_{q_{z},\omega}^{\prime*},$
(16)

where
$\bm{u}_{q_{z},\omega}^{\prime}=\left(u^{\prime}/2,-iu^{\prime}/2\right)$ and
the MA susceptibility reads

$\begin{split}\chi_{\bm{k}}=&\frac{2e^{2}\mathrm{sgn}(m)g\omega(\lambda+\mu)}{h(\lambda+2\mu)}\\\
&\times\frac{k_{y}\left(3k_{x}^{2}-k_{y}^{2}\right)\left[k_{x}^{2}f_{\bm{k},yy}-2k_{x}k_{y}f_{\bm{k},xy}+k_{y}^{2}f_{\bm{k},xx}\right]}{\left(k_{x}^{2}+k_{y}^{2}\right)^{2}}.\end{split}$
(17)

Here
$\bm{u}^{\prime}(\bm{r})=\sum_{\bm{k}}\bm{u}_{\bm{k},\omega}^{\prime}e^{i(\bm{k}\cdot\bm{r}-\omega
t)}$, $d(\bm{r})=\sum_{\bm{k}}d_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}$, and
$\partial_{i}d\partial_{j}d=\sum_{\bm{k}}f_{\bm{k},ij}e^{i\bm{k}\cdot\bm{r}}$
with
$f_{\bm{k},ij}=-\sum_{\bm{k}^{\prime}}k_{i}^{\prime}\left(k_{j}-k_{j}^{\prime}\right)d_{\bm{k}^{\prime}}d_{\bm{k}-\bm{k}^{\prime}}$,
see SM [42] for details. As demonstrated in Eq. (17), the MA susceptibility is
proportional to the acoustic frequency as well as the quantized valley Chern
number and is independent of the quasiparticle lifetime, reflecting its
topological origin. This behavior is in contrast to the conventional inverse
Faraday effect [44, 45, 46], where the magneto-optic susceptibility is
inversely proportional to the optical frequency in the dirty limit [47].
Meanwhile, the magnetization is proportional to the strain-induced charge
density and hence pseudo-magnetic field in Eq. (4), which indicates a
conversion of the emergent field into physical magnetization.

Figure 2: (a) Schematic setup for the IAFE. A periodically rippled 2D Dirac
material is placed on a piezoelectric substrate in which a circularly
polarized acoustic wave propagates in the out-of-plane direction. (b) and (c)
The renormalized MA susceptibilities in momentum and real space, respectively.
The blue and orange circles in (b) enclose the peaks at $\bm{Q}_{1}$ and
$\bm{Q}_{1}-\bm{Q}_{3}$. (d) The polar plot of the MA susceptibility peaks
$\left|\chi_{\bm{Q}_{1}}\right|$ and
$\left|\chi_{\bm{Q}_{1}-\bm{Q}_{3}}\right|$.

To be concrete, we consider a periodically rippled Dirac material with the
topography depicted by $d(\bm{r})=d_{0}\sum_{j=1}^{3}\cos\left(\bm{Q}_{j}\cdot
r\right)$ [48, 49] where $\bm{Q}_{j}=R_{3}^{j-1}\bm{Q}$ with
$\bm{Q}=Q(\cos\theta,\sin\theta)$ and the three-fold rotation operator
$R_{3}$. For $\theta=20^{\circ}$, the renormalized MA susceptibility
$\chi_{\bm{k}}/\chi_{0}$ is shown in Fig. 2(b), where the factor

$\chi_{0}=\frac{2e^{2}\mathrm{sgn}(m)g\omega(\lambda+\mu)Q^{3}d_{0}^{2}}{h(\lambda+2\mu)},$
(18)

characterizes the strength of the MA susceptibility. Interestingly, there are
twelve peaks in Fig. 2(b), including six at the wave vectors $\pm\bm{Q}_{j}$
of $d(\bm{r})$ and the others at $\pm(\bm{Q}_{j+1}-\bm{Q}_{j})$ where
$\bm{Q}_{4}\equiv\bm{Q}_{1}$. It indicates that
$\chi(\bm{r})=\sum_{\bm{k}}\chi_{\bm{k}}e^{i\bm{k}\cdot\bm{r}}$ contains
components whose period are smaller than that of $d(\bm{r})$, resulting from
the nonlinear terms in Eq. (17). In Fig. 2(c), we show
$\mathrm{Im}\chi(\bm{r})/\chi_{0}$ that exhibits the antiferromagnetic
signature. Note that $\chi(\bm{r})$ is purely imaginary such that the
magnetization
$\bm{M}(\bm{r})=\chi(\bm{r})\bm{u}_{q_{z},\omega}^{\prime}\times\bm{u}_{q_{z},\omega}^{\prime*}$
is real. Moreover, the ratio between $|\chi_{\bm{Q}_{1}}|$ and
$|\chi_{\bm{Q}_{1}-\bm{Q}_{3}}|$ can be controlled by the polar angle
$\theta$, as shown in Fig. 2(d). In particular for $\bm{Q}_{j}$ along the
zigzag (armchair) direction, $|\chi_{\bm{Q}_{1}}|$
($|\chi_{\bm{Q}_{1}-\bm{Q}_{3}}|$) vanishes, leading to the magnetization with
different periods, see SM [42] for details. For $\omega\sim 2\pi\times 1$ GHz,
$Q\sim 1$ $\mathrm{nm^{-1}}$, $d_{0}\sim 1$ nm, $\lambda=3.5$
$\mathrm{eV}/\textup{\AA}^{2}$, and $\mu=7.8$ $\mathrm{eV}/\textup{\AA}^{2}$
of monolayer hBN [50], we estimate $\chi_{0}\sim 0.4$
$\mu_{N}/\mathrm{nm^{4}}$ where $\mu_{N}$ is the nuclear magneton. For
$u^{\prime}\sim 10$ nm, it corresponds to a magnetic momentum $m_{z}\sim 1.2$
$\mu_{N}$ per unit cell, comparable to the orbital magnetic moments of phonons
in ionic materials [51, 52].

The ferromagnetism can also be realized by the IAFE, which requires the charge
density induced by pseudo-magnetic field to be uniform in space. Here we
consider a triaxial stretch of the honeycomb lattice [18, 53, 54] with
$(u_{r},u_{\theta})=u_{0}^{-1}r^{2}(\cos(3\theta-3\phi),\sin(3\theta+3\phi))$
in the polar coordinate, which leads to a uniform
$\rho=-16e^{2}h^{-1}\mathrm{sgn(m)}gu_{0}^{-1}\sin 3\phi$ in the bulk. Here
$\phi$ is the polar angle of the stretch direction. The static magnetization
generated by the circularly polarized acoustic wave is
$\bm{M}=\frac{1}{2}\rho\bm{u}^{\prime}\times\dot{\bm{u}}^{\prime}=-\frac{1}{2}\omega\rho
u^{\prime 2}\hat{\bm{z}}.$ For $\theta=90^{\circ}$and $u_{0}\sim 100$ nm, the
estimated magnetization is $M\sim 3$ $\mu_{N}/\mathrm{nm^{2}}$ for the same
circularly polarized acoustic wave above. Note that the whole system remains
charge neutral at half filling in the gapped phase. Therefore, the net charge
in the bulk must be canceled by the opposite charge at boundaries, see SM [42]
for details.

Discussion and conclusion.—Here we compare the proposed QNAHE and IAFE with
the existing nonlinear Hall effect and inverse Faraday effect. Despite the
fundamental difference between the acoustic wave and electromagnetic field as
driving forces of both effects, we further distinguish them from multiple
perspectives. So far, there are several different mechanisms of the nonlinear
Hall effect [25] that include the Berry curvature dipole [55], quantum metric
[56], disorder scattering [57], and virtual interband transition [58]. Among
them, the first three depend on the geometric and topological properties of
the electron wave function at the Fermi energy and apply only to the metal,
while the last one works for the magnetic insulator without TRS and inversion
symmetry. In contrast, QNAHE exists in the Dirac insulator with TRS but with
no inversion symmetry. The disorder scattering and virtual interband
transition, which can be substantially suppressed by the large band gap, are
irrelevant in QNAHE. Compared to the quantum valley Hall effect of the Dirac
insulator [59, 60, 61], the anomalous Hall currents from opposite valleys add
up rather than counteract in QNAHE. Therefore, it is robust against the
intervalley scattering and can be used to probe the quantized valley Chern
number. As for the inverse Faraday effect of Dirac materials, it has been
proposed through the second-order response that contains both intraband and
interband contributions in doped semimetals [62, 63, 14]. In our case, there
is no free charge carrier in the Dirac insulator and the acoustic frequency is
much smaller than the band gap. Therefore, the proposed IAFE depends only on
the intrinsic band topology characterized by the quantized valley Chern number
in the presence of TRS. We emphasize that both QNAHE and IAFE in the Dirac
insulator are independent of the quasiparticle lifetime due to their
topological origins.

In summary, we investigate the quantized nonlinear topological AE response of
the Dirac insulator that gives rise to QNAHE and IAFE. The AE conductivity and
MA susceptibility are both proportional to the quantized valley Chern number
and independent of the quasiparticle lifetime. The rectified topological
current and magnetization generated by different acoustic waves are estimated
using the realistic parameters of monolayer hBN, which are measurable within
the current experimental resolution. Owing to the large band gap of various
Dirac insulators, the proposed QNAHE and IAFE are robust and even accessible
at room temperature.

###### Acknowledgements.

Acknowledgments.—We thank Yu Zhang, Di Xiao, Kai Sun, Marcel Franz, Avadh
Saxena, Pavlo Suckhachov, Matthias Geilhufe, and Bin Yan for stimulating
discussions. The work at LANL was carried out under the auspices of the U.S.
DOE NNSA under contract No. 89233218CNA000001 through the LDRD Program, and
was supported by the Center for Nonlinear Studies at LANL, and was performed,
in part, at the Center for Integrated Nanotechnologies, an Office of Science
User Facility operated for the U.S. DOE Office of Science, under user
proposals #2018$B$$U$0010 and #2018$B$$U$0083\. AVB is supported by Quantum
CT, European Research Council under the European Union Seventh Framework
ERS-2018-SYG 810451 HERO and Knut and Alice Wallenberg Foundation KAW
2019.0068.

## References

  * Vafek and Vishwanath [2014] O. Vafek and A. Vishwanath, Annual Review of Condensed Matter Physics 5, 83 (2014).
  * T.O. Wehling and Balatsky [2014] A. B.-S. T.O. Wehling and A. Balatsky, Advances in Physics 63, 1 (2014).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Rycerz  _et al._ [2007] A. Rycerz, J. Tworzydło, and C. Beenakker, Nat. Phys. 3, 172 (2007).
  * Xiao  _et al._ [2007] D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. 99, 236809 (2007).
  * Schaibley  _et al._ [2016] J. R. Schaibley, H. Yu, G. Clark, P. Rivera, J. S. Ross, K. L. Seyler, W. Yao, and X. Xu, Nat. Rev. Mater. 1, 16055 (2016).
  * Vitale  _et al._ [2018] S. A. Vitale, D. Nezich, J. O. Varghese, P. Kim, N. Gedik, P. Jarillo-Herrero, D. Xiao, and M. Rothschild, Small 14, 1801483 (2018).
  * Kalameitsev  _et al._ [2019] A. V. Kalameitsev, V. M. Kovalev, and I. G. Savenko, Phys. Rev. Lett. 122, 256801 (2019).
  * Sonowal  _et al._ [2020] K. Sonowal, A. V. Kalameitsev, V. M. Kovalev, and I. G. Savenko, Phys. Rev. B 102, 235405 (2020).
  * Wan  _et al._ [2024] J.-L. Wan, Y.-L. Wu, K.-Q. Chen, and X.-Q. Yu, Phys. Rev. B 109, L161101 (2024).
  * Sukhachov and Rostami [2020] P. O. Sukhachov and H. Rostami, Phys. Rev. Lett. 124, 126602 (2020).
  * Zhao  _et al._ [2022] P. Zhao, C. H. Sharma, R. Liang, C. Glasenapp, L. Mourokh, V. M. Kovalev, P. Huber, M. Prada, L. Tiemann, and R. H. Blick, Phys. Rev. Lett. 128, 256601 (2022).
  * Bhalla  _et al._ [2022] P. Bhalla, G. Vignale, and H. Rostami, Phys. Rev. B 105, 125407 (2022).
  * Liang  _et al._ [2021] L. Liang, P. O. Sukhachov, and A. V. Balatsky, Phys. Rev. Lett. 126, 247202 (2021).
  * Suzuura and Ando [2002] H. Suzuura and T. Ando, Phys. Rev. B 65, 235412 (2002).
  * Mañes [2007] J. L. Mañes, Phys. Rev. B 76, 045430 (2007).
  * Vozmediano  _et al._ [2010] M. A. Vozmediano, M. Katsnelson, and F. Guinea, Physics Reports 496, 109 (2010).
  * Guinea  _et al._ [2010] F. Guinea, M. I. Katsnelson, and A. Geim, Nature Physics 6, 30 (2010).
  * Levy  _et al._ [2010] N. Levy, S. Burke, K. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. C. Neto, and M. F. Crommie, Science 329, 544 (2010).
  * Cazalilla  _et al._ [2014] M. A. Cazalilla, H. Ochoa, and F. Guinea, Phys. Rev. Lett. 113, 077201 (2014).
  * Settnes  _et al._ [2016a] M. Settnes, S. R. Power, M. Brandbyge, and A.-P. Jauho, Phys. Rev. Lett. 117, 276801 (2016a).
  * Li  _et al._ [2020] S.-Y. Li, Y. Su, Y.-N. Ren, and L. He, Phys. Rev. Lett. 124, 106802 (2020).
  * Parmenter [1953] R. H. Parmenter, Phys. Rev. 89, 990 (1953).
  * Weinreich and White [1957] G. Weinreich and H. G. White, Phys. Rev. 106, 1104 (1957).
  * Du  _et al._ [2021] Z. Du, H.-Z. Lu, and X. Xie, Nature Reviews Physics 3, 744 (2021).
  * Tokman and Pozdnyakova [2013] I. D. Tokman and V. I. Pozdnyakova, The European Physical Journal B 86, 207 (2013).
  * Shabala and Geilhufe [2024] N. Shabala and R. M. Geilhufe, arXiv preprint arXiv:2405.09538 (2024).
  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Phys. Rev. Lett. 45, 494 (1980).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * Haldane [1988] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
  * Bernevig and Zhang [2006] B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
  * Konig  _et al._ [2007] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, _et al._ , Science 340, 167 (2013).
  * Vaezi  _et al._ [2013] A. Vaezi, N. Abedpour, R. Asgari, A. Cortijo, and M. A. H. Vozmediano, Phys. Rev. B 88, 125406 (2013).
  * Sela  _et al._ [2020] E. Sela, Y. Bloch, F. von Oppen, and M. B. Shalom, Phys. Rev. Lett. 124, 026602 (2020).
  * Droth  _et al._ [2016] M. Droth, G. Burkard, and V. M. Pereira, Phys. Rev. B 94, 075404 (2016).
  * Rostami  _et al._ [2018] H. Rostami, F. Guinea, M. Polini, and R. Roldán, npj 2D Materials and Applications 2, 15 (2018).
  * Rodriguez-Lopez and Cortijo [2018] P. Rodriguez-Lopez and A. Cortijo, Phys. Rev. B 97, 235128 (2018).
  * White and Voltmer [1965] R. M. White and F. W. Voltmer, Applied physics letters 7, 314 (1965).
  * Biryukov  _et al._ [1995] S. V. Biryukov, Y. V. Gulyaev, V. V. Krylov, and V. P. Plessky,  _Surface acoustic waves in inhomogeneous media_ , Vol. 20 (Springer, 1995).
  * [42] See Supplemental Materials for (1) symmetry properties of the acoustoelectric conductivity, (2) strain-induced charge density, (3) derivation of the magnetoacoustic susceptibility, (4) polar angle dependence of the magnetoacoustic susceptibility, which includes Refs. [64, 37, 65, 43].
  * Guinea  _et al._ [2008] F. Guinea, B. Horovitz, and P. Le Doussal, Phys. Rev. B 77, 205421 (2008).
  * Pitaevskii [1961] L. Pitaevskii, Sov. Phys. JETP 12, 1008 (1961).
  * van der Ziel  _et al._ [1965] J. P. van der Ziel, P. S. Pershan, and L. D. Malmstrom, Phys. Rev. Lett. 15, 190 (1965).
  * Pershan  _et al._ [1966] P. S. Pershan, J. P. van der Ziel, and L. D. Malmstrom, Phys. Rev. 143, 574 (1966).
  * Battiato  _et al._ [2014] M. Battiato, G. Barbalinardo, and P. M. Oppeneer, Phys. Rev. B 89, 014413 (2014).
  * Vázquez de Parga  _et al._ [2008] A. L. Vázquez de Parga, F. Calleja, B. Borca, M. C. G. Passeggi, J. J. Hinarejos, F. Guinea, and R. Miranda, Phys. Rev. Lett. 100, 056807 (2008).
  * Guinea  _et al._ [2009] F. Guinea, B. Horovitz, and P. Le Doussal, Solid State Communications 149, 1140 (2009).
  * Jung  _et al._ [2015] J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nature communications 6, 6308 (2015).
  * Juraschek  _et al._ [2017] D. M. Juraschek, M. Fechner, A. V. Balatsky, and N. A. Spaldin, Phys. Rev. Mater. 1, 014401 (2017).
  * Juraschek and Spaldin [2019] D. M. Juraschek and N. A. Spaldin, Phys. Rev. Mater. 3, 064405 (2019).
  * Neek-Amal  _et al._ [2013] M. Neek-Amal, L. Covaci, K. Shakouri, and F. M. Peeters, Phys. Rev. B 88, 115428 (2013).
  * Settnes  _et al._ [2016b] M. Settnes, S. R. Power, and A.-P. Jauho, Phys. Rev. B 93, 035456 (2016b).
  * Sodemann and Fu [2015] I. Sodemann and L. Fu, Phys. Rev. Lett. 115, 216806 (2015).
  * Gao  _et al._ [2014] Y. Gao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 112, 166601 (2014).
  * Du  _et al._ [2019] Z. Du, C. Wang, S. Li, H.-Z. Lu, and X. Xie, Nature Communications 10, 3047 (2019).
  * Kaplan  _et al._ [2023] D. Kaplan, T. Holder, and B. Yan, Nature communications 14, 3053 (2023).
  * Yao  _et al._ [2009] W. Yao, S. A. Yang, and Q. Niu, Phys. Rev. Lett. 102, 096801 (2009).
  * Gorbachev  _et al._ [2014] R. Gorbachev, J. Song, G. Yu, A. Kretinin, F. Withers, Y. Cao, A. Mishchenko, I. Grigorieva, K. S. Novoselov, L. Levitov, _et al._ , Science 346, 448 (2014).
  * Shimazaki  _et al._ [2015] Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watanabe, T. Taniguchi, and S. Tarucha, Nature Physics 11, 1032 (2015).
  * Tokman  _et al._ [2020] I. D. Tokman, Q. Chen, I. A. Shereshevsky, V. I. Pozdnyakova, I. Oladyshkin, M. Tokman, and A. Belyanin, Phys. Rev. B 101, 174429 (2020).
  * Gao  _et al._ [2020] Y. Gao, C. Wang, and D. Xiao, arXiv preprint arXiv:2009.13392 (2020).
  * Ando [2006] T. Ando, Journal of the physical society of Japan 75, 124701 (2006).
  * Bernevig and Hughes [2013] B. A. Bernevig and T. L. Hughes,  _Topological insulators and topological superconductors_ (Princeton university press, 2013).

in 1,…,0 See pages \x, of supplement.pdf