# Bilayer Hubbard model: Analysis based on the fermionic sign problem

Yingping Mou  Beijing Computational Science Research Center, Beijing 100193,
China  R. Mondaini<EMAIL_ADDRESS>Beijing Computational Science
Research Center, Beijing 100193, China  R.T. Scalettar
<EMAIL_ADDRESS>Department of Physics, University of California,
Davis, California 95616, USA

###### Abstract

The bilayer Hubbard model describes the antiferromagnet to spin singlet
transition and, potentially, aspects of the physics of unconventional
superconductors. Despite these important applications, significant aspects of
its ‘phase diagram’ in the interplane hopping $t_{\perp}$ versus on-site
interaction $U$ parameter space, at half filling, are largely in disagreement.
Here we provide an analysis making use of the average sign of weights over the
course of the importance sampling in quantum Monte Carlo simulations to
resolve several central open questions. Specifically, this metric of the
weights clarifies the finite-sized metallic regimes at small $U$. Furthermore,
at strong interactions, it points to the existence of a crossover from a
correlated to uncorrelated band insulator not yet explored in a variety of
existing, unbiased numerical methods. Our paper demonstrates the versatility
of using properties of the weights in quantum Monte Carlo simulations to
reveal important physical characteristics of the models under study.

##  I Introduction

In the wake of the discovery of high-temperature superconductors, and the
unlikeliness of their description within a conventional electron-phonon
mechanism [1], strongly interacting models that could explain the physical
mechanisms occurring in cuprates had a surge of investigation. Among these,
multilayer geometries such as the bilayer Heisenberg model [2, 3], are
essential to understand the robust (i.e., finite temperature)
antiferromagnetic ordering observed in undoped materials. The bilayer $t$-$J$
[4, 5] and Hubbard models [6, 7, 8] allowed the study of the interplay of
itinerant electrons and (short-ranged) magnetic ordering in the presence of
hole doping, and hence spin fluctuation mediated pairing. In the latter,
original studies have pointed out the possibility of a nodeless $d$-wave
pairing, where the gap has opposite signs in the bonding and anti-bonding
Fermi surfaces, and that interplane hybridization weakens in-plane
superconducting correlations.

Due to the presence of the sign problem [9, 10, 11] in the doped regime [12],
investigations using quantum Monte Carlo (QMC) simulations had most success
studying the half-filled case [13, 14], which allows the understanding of how
global long-range magnetic ordering takes place at sufficiently small
interplane hybridizations. More recently, large scale ground-state QMC
calculations [15] have clarified the absence of metallicity at finite values
of the interactions, as initially suggested to occur [16, 17, 13, 18, 19], and
further corroborated the existence of a magnetic transition in the 3D
Heisenberg universality class as the interplane hybridization was increased,
similar to that of the bilayer Heisenberg model [20]. Dynamical cluster
approximation calculations have also examined the possibility of the
enhancement of superconductivity as one of the bands approaches the Lifshitz
transition, and its implications for heavily electron-doped FeSe-derived
superconductors [21].

In this paper, we revisit the phase diagram of the half-filled bilayer Hubbard
model using the finite-temperature determinant quantum Monte Carlo (DQMC)
method [22, 23], with a goal of establishing its different phase boundaries in
a way that takes the average sign of partial weights in the sampling as a
minimal correlator. This metric, recently used to understand criticality in a
variety of quantum models [24], further allows one to unveil a subtle and
often unappreciated crossover from a correlated to uncorrelated band-
insulating regime at large interplane hybridizations.

Figure 1: (a) Schematic illustration of the Hamiltonian, with relevant
parameters annotated. (b) The contour plot of the average spin resolved sign
$\langle{\cal S}_{\sigma}\rangle$ in the space of parameters $(t_{\perp},U)$;
(c) the equivalent for the antiferromagnetic structure factor $S_{\rm AFM}$.
All data are extracted at temperatures $T/t=1/20$ in a $12\times 12$ bilayer.
Markers along the $U/t=0$ axis in (b) and (c) depict the values of
$t_{\perp}/t$ that result in a nesting condition with a wavevector $(\pi,\pi)$
for this lattice size. Examples of such nestings are shown in (d) and (e) for
$t_{\perp}/t=1.73$ and 3.0, respectively.Lines in (b) and (c) depict the
results of Ref. 15 marking the magnetic transition obtained via the scaling of
the antiferromagnetic order parameter at $T=0$.

##  II Model

The Hamiltonian of an $L\times L$ bilayer reads

$\displaystyle\hat{H}=$ $\displaystyle-$ $\displaystyle t\sum_{\langle
ij\rangle\ell\,\sigma}(\hat{c}^{\dagger}_{i\ell\sigma}\hat{c}^{\phantom{\dagger}}_{j\ell\sigma}+{\rm
H.c.})-t_{\perp}\sum_{i\,\sigma}(\hat{c}^{\dagger}_{i0\sigma}\hat{c}^{\phantom{\dagger}}_{i1\sigma}+{\rm
H.c.})$ (1) $\displaystyle+$ $\displaystyle
U\sum_{i\ell}\left(\hat{n}_{i\ell\uparrow}-\frac{1}{2}\right)\left(\hat{n}_{i\ell\downarrow}-\frac{1}{2}\right)-\mu\sum_{i,l,\sigma}\hat{n}_{il\sigma},$

where $\hat{c}^{\dagger}_{i\ell\sigma}$ is the fermionic creation operator on
site $i$ of the layer $\ell$ ($\ell=0,1$) with spin $\sigma$
($\sigma=\uparrow,\downarrow$), and $\hat{n}_{i\ell\sigma}$ is the
corresponding number density operator; $t$ and $t_{\perp}$ quantify the
nearest-neighbor intra- and interplane hopping amplitudes, while $U$ is the
strength of the local repulsive interactions with chemical potential $\mu$
controlling the electronic density [Fig. 1(a)].

In the non-interacting limit ($U=0$) at half-filling ($\mu=0$), the system
undergoes a metal-to-band insulator transition as $t_{\perp}/t>4$, with a gap
opening between the bonding ($-$) and anti-bonding ($+$) bands,
$\varepsilon_{0}^{\pm}({\bf k})=-2t[\cos(k_{x})+\cos(k_{y})]\pm t_{\perp}$,
whose size is $t_{\perp}-4t$. In the opposite, strongly interacting, case
($U\gg t,t_{\perp}$), the Hamiltonian at $\mu=0$ is equivalent to the bilayer
Heisenberg model [2, 20], which displays a magnetic transition from an
antiferromagnetic ordered bilayer to a quantum disordered phase featuring
interplane singlets. Mapping the spin exchange interactions $J$ to the hopping
energy scales at this limit, gives a critical hybridization
$t_{\perp}^{c}/t=1.588$ separating these two regimes.

For the generic $U\neq 0$ case, we solve Eq.(1) by making use of DQMC, in
which the introduction of a real-space imaginary-time auxiliary field
$\\{s_{i\tau}\\}$ decouples the interactions, allowing the fermionic
integration to be taken exactly. As a result, the partition function is
written in terms of the product of weights for each fermionic flavor $\sigma$,
${\cal Z}=\sum_{s}\prod_{\sigma}w_{\sigma}(\\{s_{i\tau}\\})$, where the field
$\\{s_{i\tau}\\}$ is being summed. Instead of solving for all configurations
$\\{s_{i\tau}\\}$, importance sampling is performed while observing the
convergence of physical observables. The single, controllable approximation
used is the imaginary-time discretization $\Delta\tau$ which we take as 0.1
throughout.

The form of the partition function reveals a peculiarity of the method: the
weights being summed are not positive-definite. In fact, in a large class of
problems of interest this ‘sign problem’ precludes the accurate computation of
physical quantities in the most interesting parts of their phase diagrams [25,
9]. At half-filling, however, due to the bipartite structure of the lattice,
there is no sign problem [23]. That is, the product of the signs of the
weights in the Monte Carlo sampling is positive regardless of the
configuration $\\{s_{i\tau}\\}$ of the Hubbard-Stratonovich (HS) field.
However, even though the total sign of the weights is always positive, the
sign of individual ones are not. Indeed, recent results have demonstrated that
the average sign of individual weights, ${\langle\cal
S}_{\sigma}\rangle\equiv\langle{\rm sgn}(w_{\sigma})\rangle$, are directly
related to the physics of the Hamiltonian under investigation, following
scaling laws similar to those for physical observables in the vicinity of
quantum phase transitions [24].

Figure 2: (a)-(d) The single-particle excitation spectrum $A_{\bf k}(\omega)$
across a path connecting high-symmetry points ($\Gamma-X-M-\Gamma$) in the
Brillouin zone, for the indicated hybridization amplitudes In (a),
$t_{\perp}/t=1.5$, and (d), $t_{\perp}/t=4.0$, the single-particle gap
$\Delta_{\rm sp}$ is marked, exhibiting, respectively, a direct and indirect
behavior, respectively. (f)-(h) Dynamical spin structure factor $S({\bf
k},\omega)$ at different $t_{\perp}/t$ as marked, with a momentum integrated
version shown in (e). Data are obtained at temperatures $T/t=1/20$ in an
$L=12$ bilayer with $U/t=10$.

##  III The sign problem phase diagram

We start by studying the sign phase diagram of the bilayer Hubbard model in
Fig. 1(b), focusing on the regime $U\lesssim t$. Lobes where $\langle{\cal
S}_{\sigma}\rangle$ is close to 1 appear at small interaction strengths, while
in a broad region of parameters $\langle{\cal S}_{\sigma}\rangle\to 0$. Not
coincidentally, the boundaries of such lobes in the $U/t\to 0$ limit match the
loci where a perfect nesting of the Fermi surfaces, $\varepsilon_{0}^{+}({\bf
k}+{\bf Q})=-\varepsilon_{0}^{-}({\bf k})$, occur in a finite lattice [see
Figs. 1(d) and 1(e) for two examples] when spanning the interplane hybrization
$t_{\perp}/t$ [15]. These lobes are finite-size effects which vanish in the
thermodynamic limit, where the spacing between nesting conditions similarly
vanishes other system sizes are shown in Appendix A. Moreover, the ${\bf
Q}=(\pi,\pi)$ nesting suggests that the system is unstable toward
antiferromagnetic order, and ensuing Mott insulating behavior is the most
probable scenario, as further pointed out by mean-field calculations [15]. The
immediate conclusion is that the regions at small $U/t$ with $\langle{\cal
S}_{\sigma}\rangle\simeq 1$ identify the metallic regime observed in finite
lattices over a variety of numerical methods [16, 17, 13, 18, 19].

Going beyond small interactions, studies that use $\langle{\cal
S}_{\sigma}\rangle$ to track the phase transitions [24] have shown that the
region with a vanishing spin-resolved sign can be identified, in related
models and when approaching the thermodynamic limit, with either the
magnetically ordered regime or the one yielding a Mott insulator. This is the
case for the SU(2) Hubbard model in the honeycomb lattice, for example, where
both phases are known to concomitantly occur [26, 27, 28]. Here, in the SU(2)
bilayer Hubbard model, however, by overlaying in Fig. 1(b) the accurate
numerical results from Ref. [15] for the critical interaction strength
$U_{c}(t_{\perp})$ that leads to the onset of magnetism in the $T\to 0$ limit,
we demonstrate that $\langle{\cal S}_{\sigma}\rangle\to 0$ is not tracking the
magnetic phase transition. Our results resolve the transition (or crossover)
it probes, and, in doing so, provide insight into the transition between
uncorrelated and correlated band insulating regimes.

This becomes more evident if contrasting Figs. 1(b) and 1(c). By plotting the
antiferromagnetic structure factor $S_{\rm
AFM}\equiv(1/2L^{2})\sum_{i,j}(-1)^{i+j}\langle(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow})(\hat{n}_{j\uparrow}-\hat{n}_{j\downarrow})\rangle$,
in a relatively large lattice ($L=12$) [Fig. 1(c)], the critical points
$U_{c}(t_{\perp})$ systematically border the regime where $S_{\rm AFM}$ is
large. The average spin-resolved sign, on the other hand, departs from zero at
larger interplane hybridizations [Fig. 1(b)]. In particular, that the magnetic
transition is not accompanied by a Mott transition in the bilayer Hubbard
model was initially shown within cluster dynamical mean-field (DMFT)
calculations [17], which demonstrated the existence of a paramagnetic Mott
insulator preceding the onset of a band insulating state at larger $t_{\perp}$
111Exact diagonalization results in Appendix G provide further evidence that
the magnetic and the Mott insulating transitions are not aligned at large
interaction strengths in the bilayer Hubbard model..

##  IV The Mott insulator-Band insulator crossover

Differentiating between Mott and band insulators at sufficiently large
$t_{\perp}$ is challenging. This difficulty has been illustrated not only in
theoretical studies of model Hamiltonians, but also in the experimental
characterization of certain transition-metal dichalcogenides [30], which
exhibit a competition of on-site Coulomb repulsion and interlayer hopping in
its layered structure. As both phases naturally manifest a finite gap
$\Delta_{\rm sp}$ for single-particle excitations, either separating the upper
and lower Hubbard bands in the Mott phase or the bonding and anti-bonding
bands in the band insulating regime, a useful distinguishing characteristic is
provided by the trend of $\Delta_{\rm sp}$ with growing interplane
hybridization [17] as well as the direct or indirect nature (in momentum) of
$\Delta_{\rm sp}$.

To implement this approach, we start by reporting the single-particle spectral
function $A_{\bf k}(\omega)$ in Fig. 2, calculated via the stochastic analytic
continuation of imaginary-time dependent QMC data [31],

${\bm{G}}({\bf
k},\tau)=\int\frac{d\omega}{\pi}\frac{e^{-\omega\tau}}{1+e^{-\beta\omega}}A_{\bf
k}(\omega),$ (2)

where ${\bm{G}}({\bf k},\tau)$ is the space Fourier transform of the imaginary
time ($\tau$) displaced Green’s function ${\bm{G}}({\bf
k},\tau)=\langle\Psi_{{\bf r+R},\sigma}(\tau)\Psi^{{\dagger}}_{{\bf
r},\sigma}(0)\rangle$, where the creation operator is given as
$\Psi^{{\dagger}}_{{\bf r},\sigma}=({\hat{c}}^{{\dagger}}_{{\bf
r}0,\sigma},{\hat{c}}^{{\dagger}}_{{\bf r}1,\sigma})$ with the creation
operator of the electron at unit cell ${\bf r}$ and sublattice $0,1$ with spin
$\sigma$.

By focusing on large interactions, $U/t=10$, we notice that at small
$t_{\perp}$ [$t_{\perp}\lesssim 2t$, see Fig. 2(a)], the single-particle
excitation bands display a direct gap at the $X=(\pi,0)$ point of the
Brillouin zone as in the single-layer Hubbard model, whereas for larger
interlayer hoppings ($t_{\perp}\gtrsim 2t$), the gap becomes indirect,
connecting the $\Gamma$-$M$ [$(0,0)$-$(\pi,\pi)$] points instead, Figs. 2(c)
and 2(d). Moreover, for values of $t_{\perp}\approx 3t$, the indirect gap
reaches a minimum, whereupon increasing hybridization leads to a larger
$\Delta_{\rm sp}$. This is summarized in Fig. 3(b) for different system sizes
additional details in the $A_{\bf k}(\omega)$ results are described in
Appendix B, and an independent scheme (with similar results) that bypasses the
analytic continuation is given in Appendix C.

Figure 3: (a) A cut along $U/t=10$ of the spin-resolved average sign
$\langle{\cal S}_{\sigma}\rangle$ with growing interplane hopping $t_{\perp}$
and various system sizes. The inset highlights the significant drift in the
upturn location $t_{\perp}/t$ of $\langle{\cal S}_{\sigma}\rangle$ as one
decreases the temperature $T$ (or increases the inverse temperature
$\beta=1/T$) for an $L=12$ bilayer. (b) The corresponding single-particle gap
$\Delta_{\rm sp}$ extracted by the minimum gap of single-particle excitations
across the whole Brillouin zone, at comparable system sizes; also shown the
gap to spin excitations $\Delta_{\rm s}$, finite once the AFM order is absent.
This is also demonstrated by the vanishing of the normalized (by the system
size) AFM structure factor, shown in (c), in the thermodynamic limit at
$t_{\perp}=2.05(5)t$ see first vertical shaded region in (b). The second
shaded region at $t_{\perp}=5.75(25)t$ describes the onset of $\Delta_{\rm
sp}\simeq\Delta_{s}$, also seen in the zoom-in of the saturation of
$\langle{\cal S}_{\sigma}\rangle$ in (d). Data are extracted at temperatures
$T/t=1/20$.

This analysis allows us to divide the dependence of $\Delta_{\rm sp}$ on the
hybridization at fixed interaction strength in three classifications: (i) a
gap that is largely robust to the increase of $t_{\perp}$ up until
$t_{\perp}/t\simeq 1.5-2$, (ii) a decrease of the gap with $t_{\perp}$ as the
hybridization approaches $t_{\perp}/t\simeq 3$, and (iii) a growth of
$\Delta_{\rm sp}$ for values $t_{\perp}/t\gtrsim 3$. The observed trend of the
gap $d\Delta_{\rm sp}/dt_{\perp}<0$ ($d\Delta_{\rm sp}/dt_{\perp}>0$) was used
before [17] in cluster-DMFT results to discern the Mott (band-insulating)
phases. Consequently, the three regimes we describe can be interpreted,
respectively, as the antiferromagnetic Mott insulator, paramagnetic Mott
insulator 222We are currently unaware of other possible orders that may
manifest in the described paramagnetic Mott insulator. A rung bond-ordered
wave is not a good candidate – see Appendix D., and band insulator.
Noticeably, quantitatively similar conclusions about the gap change and the
Mott-to-band insulator transition were obtained in smaller clusters using
dynamical cluster approximation [33], albeit at more modest interactions
($U/t=6$).

##  V Characterizing band insulators

The classification of the onset of a band insulator via the evolution of
$\Delta_{\rm sp}$ with $t_{\perp}$ is however incomplete. It has been argued
that a true (or uncorrelated) band insulator is characterized by identical
gaps for particle and spin excitations [34]. To verify this, we compute the
dynamical spin structure factor, $S(\textbf{k},\omega)$ via the inversion of
the integral equation,

$\langle\hat{S}({\bf
k},\tau)\hat{S}({\bf-k},0)\rangle=\int\frac{d\omega}{\pi}\frac{e^{-\tau\omega}}{1-e^{-\beta\omega}}\chi^{\prime\prime}({\bf
k},\omega),$ (3)

where $\chi^{\prime\prime}({\bf k},\omega)$ is the dynamical spin
susceptibility, which allows the extraction of $S({\bf
k},\omega)=\chi^{\prime\prime}({\bf k},\omega)/(1-e^{-\beta\omega})$. Figures
2(f)-2(h) reports the dynamical spin structure factor for a representative
path in the Brillouin zone. A finite spin gap $\Delta_{s}$ appears at values
$t_{\perp}/t\gtrsim 2$ at the $M=(\pi,\pi)$ point, in agreement with the
disappearance of antiferromagnetic long-range order, as also seen in the $1/L$
extrapolation of the normalized equal-time antiferromagnetic spin structure
factor, $S_{\rm AFM}=\frac{1}{2L^{2}}\sum_{i,j}{\rm exp}\\{{\rm
i}(\pi,\pi)\cdot({\bf r}_{i}-{\bf
r}_{j})\\}\langle(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})(\hat{n}_{j,\uparrow}-\hat{n}_{j,\downarrow})\rangle$
[Fig. 3(c)]. Larger interplane hybridizations lead to a much reduced momentum
dependence in $S({\bf k},\omega)$; a momentum-integrated version is shown in
Fig. 2(e), and a compilation of the spin gap extracted from $S(\omega)$ is
presented in Fig. 3(b). Both gaps, $\Delta_{\rm sp}$ and $\Delta_{s}$, acquire
similar values at $t_{\perp}/t\approx 5.8$, thus marking the onset of the
uncorrelated band-insulating phase.

We are now in position to analyze how $\langle{\cal S}_{\sigma}\rangle$
captures different crossovers, as reported in Fig. 3(a). Despite not
insignificant system-size dependence, these data indicate that the upturn of
the average spin-resolved sign at sufficiently low temperatures occurs close
to the Mott insulator-band insulator crossover at around $t_{\perp}\approx 3t$
[as originally seen in Fig. 1(b)]. In turn, $\langle{\cal S}_{\sigma}\rangle$
saturates at one precisely when the uncorrelated band-insulating regime takes
place. Figure 3(d) shows a detailed zoom-in of this $\langle{\cal
S}_{\sigma}\rangle\to 1$ approach. Such saturation is also seen in quantifiers
of the typically sampled fields in DQMC, e.g., the average Hamming distance
[35] –see Appendix F. Although a scaling analysis based on the average sign of
the weights $\langle{\cal S}_{\sigma}\rangle$ [24] is elusive, likely due to
the absence of an intrinsic order parameter characterizing either phase, other
models where a band-insulating phase takes place, as the square lattice ionic
Hubbard model [36, 37, 38, 39, 40], also exhibit a convergence of
$\langle{\cal S}_{\sigma}\rangle$ towards one, as we similarly observe here
[11, 24]. Further characterization of this crossover is shown in Appendix E
which describes several local spin correlators.

##  VI Summary and outlook

By using a combination of quasiparticle excitation gaps with the average sign
of the weight of one fermionic flavor $\langle{\cal S}_{\sigma}\rangle$ in the
importance sampling of DQMC simulations, we identified the different phases of
the bilayer Hubbard model at half filling. Beyond the clarification of the
system size-influenced metallic regimes at small interactions, one of our main
results using this tracker is the onset of an uncorrelated band-insulating
regime at large interplane hybridizations. While this transition is likely a
crossover, i.e., without an associated order parameter, the saturation of
$\langle{\cal S}_{\sigma}\rangle$ at one coincides with the regime where spin
and single-particle excitations are comparable. As a result, in both
uncorrelated phases at finite lattices, metallic and band insulating,
$\langle{\cal S}_{\sigma}\rangle\to 1$. That $\langle{\cal S}\rangle$ might
signal transitions in model Hamiltonians is consistent with its role as a
necessary ingredient to compute any physical observable [23], thus inherently
revealing details of the physics at play. Moreover, for the specific case of
the crossover between band insulating phases, it has the added advantage that
it does not require the quantification of, often expensive, time-displaced
correlation functions, neither analytic continuation of the data.

A testament of the relevance of these results is provided by experiments
involving ultracold atoms trapped in optical lattices, which have recently
succeeded in emulating the bilayer Hubbard model [41, 42], opening a further,
and highly controllable, realization for precision investigations of its
different phases.

Going beyond half filling, our exploration of the sign problem in this model
may allow one to tackle the regime which is mostly relevant to the doped
cuprates, where bilayer features have been observed [43] and numerically
interpreted [44], culminating in the investigation of the potential onset of a
finite-temperature superconducting transition governed by a Kosterlitz-
Thouless form. Such analysis has been recently carried out for the case of a
single layer [24]. Its extension to a bilayer model might reveal the impact of
the interplane hybridization on the finite critical temperature $T_{c}$, as
well as on the dominant pairing channel before this transition takes place
[45, 46, 47, 48].

##  VII Acknowledgements

Y.M. was supported by Grant No. 12047507 funded by the Natural Science
Foundation of China (NSFC). R.M. acknowledges discussions with H. Guo, and
support from the NSFC Grants No. U1930402, No. 12050410263, No. 12111530010,
No. 11974039, and No. 12222401. R.T.S. was supported by Grant DE‐SC0014671
funded by the U.S. Department of Energy, Office of Science. Illuminating
discussions with M. Randeria are gratefully appreciated. Computations were
performed with the QUEST toolbox [49] on the Tianhe-2JK at the Beijing
Computational Science Research Center.

Figure 4:  (a) [(b)] The average spin-resolved sign $\langle{\cal
S}_{\sigma}\rangle$ [derivative of the double occupancy with respect to
$t_{\perp}$, $d\langle\hat{n}_{\uparrow\downarrow}\rangle/dt_{\perp}$] on the
$T$ versus $U/t$ map for a $4\times 4$ bilayer. (c) and (d) show the same for
an $L=8$ lattice. As in Fig. 1, markers along the $U/t=0$ axis denote nesting
conditions for the corresponding lattice size, and data are extracted at
temperature $T/t=1/20$.

##  Appendix A Local quantities and different system sizes

From the discussion presented in the main text, the location of the system
size-influenced metallic regions is already evidenced by the regimes where the
$\langle{\cal S}_{\sigma}\rangle$ converges to one. To make this connection
clearer, and observe how they extend as the lattice size is reduced, we show
in Figs. 4(a) and 4(c), the ‘sign phase diagram’ for $L=4$ and 8,
respectively. The lobe boundaries in the noninteracting limit are again given
by the values of $t_{\perp}/t$ where a nesting condition occur, even more
clearly as the number of such matches is proportional to the system size. A
systematic shrinking of these regions takes place with growing $L$.

While the weight of the configurations in the sampling might initially be
thought to be an artifact of the quantum simulation algorithm, it directly
captures the behavior of physical observables, including the double occupancy,
$\langle\hat{n}_{\uparrow\downarrow}\rangle\equiv(1/2L^{2})\sum_{i}\langle\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\rangle$.
In an insulator-to-metal (metal-to-insulator) transition, the double occupancy
increases (decreases) with changes in the parameter driving the evolution, for
example, the interplane hybridization. Figures 4(b) and 4(d) show that the
derivative of the double occupancy with respect to $t_{\perp}$,
$d\langle\hat{n}_{\uparrow\downarrow}\rangle/dt_{\perp}$, reflects this
prediction, and identifies the metallic region boundaries in a finite lattice,
in direct agreement with the average sign analysis.

Figure 5: The spectral function $A_{\bf k}(\omega)$ (a)-(d) and the single-
particle density of states $N(\omega)$ (e)-(h) for various interplane hoppings
as marked. For values $t_{\perp}/t\gtrsim 2$ the gap turns indirect,
connecting the two bands at the conduction and valence bands at $(\pi,\pi)$
and $(0,0)$, respectively. The density of states at large $t_{\perp}/t$’s
displays a collection of features connected to the incoherent superposition of
interplane singlet-states. Data is obtained for $U/t=10$, $T/t=1/20$ in an
$L=12$ bilayer.

##  Appendix B QMC results: Spectral function and DOS

Following the discussion in the main text that establishes that the magnetic
transition is accompanied by a transition from a direct-to-indirect gap of the
single-particle excitations, we display in Figs. 5(a)-5(d) the momentum-
resolved $A_{\bf k}(\omega)$ for a range of $t_{\perp}$ values. To contrast
these results, we further show the full density of states,
$N(\omega)=\sum_{\bf k}A_{\bf k}(\omega)$ in Figs. 5 (e)-5(h).

For the spectral function, the change of the gap type, from direct to
indirect, is seen to happen at values $t_{\perp}/t\gtrsim 2$, with a direct
$[(\pi,0)-(\pi,0)]$ gap giving way to $[(\pi,\pi)-(0,0)]$, closely following
the magnetic transition at $U/t=10$ [15]. Deep in the nonmagnetic phase the
$(\pi,\pi)$ contribution to the excitations departs from the Fermi energy, and
thus becomes less relevant at low energies. In turn, the density of states,
which is always gapped when increasing the interplane hoppings, displays a
variety of incoherent peaks that are a characteristic of the superposition of
singlet-states formed across the bilayer in the regime $t_{\perp}/t\gg 1$. In
the main text, Fig. 3(b), the gaps were extracted via the energy difference
given by the location of the maximum values of the $\delta$-like functions
closest to the Fermi energy at $\omega=0$.

Figure 6: (a),(b) The imaginary-time ($\tau$) dependent Green’s function
$G(\bf k,\tau)$ at $U/t=10$, $t_{\perp}/t=4.5$ for $L=12$. Inset: $G({\bf
k},\tau)$ on semi-logarithmic scales. (c) The single-particle gap $\Delta_{\rm
sp}$ extracted by the minimum gap of fitting ${\rm ln}G({\bf
k},\tau)\sim\pm\omega_{\mp}({\bf k})\tau$ across the whole Brillouin zone (see
text).

##  Appendix C Directly extracting the single-particle gap

In the main text, we make use of stochastic analytic continuation to invert
the integral equation [Eq. (2)] that allows one to retrieve the spectral
function $A_{\bf k}(\omega)$ (and consequently the single-particle gap) via
the examination of the smallest gap across the allowed momentum points,
$\Delta_{\rm sp}({\bf k})$. Here, we employ a direct approach that complements
this procedure. The single-particle gap can be similarly obtained by noticing
that $G({\bf k},\tau)\propto e^{\pm\omega_{\mp}({\bf k})\tau}$, for imaginary-
time $\tau$ sufficiently away from its limits at $\tau=0$ and $\tau=\beta$. In
particular, for the case of $\tau<\beta/2$ ($\tau>\beta/2$) one obtains the
gap $\omega_{+}$ ($\omega_{-}$) above (below) the Fermi energy at $\mu=0$. As
a result, the single-particle gap is extracted as $\Delta_{\rm sp}=\min_{\bf
k}[\omega_{+}({\bf k})]+\min_{\bf k}[\omega_{-}({\bf k})]$. Figures 6(a) and
6(b) exemplify the $\tau$ dependence of the imaginary-time displaced Green’s
functions at two momentum points, ${\bf k}=(0,0)$ and $(\pi,\pi)$,
respectively. A least-squares fitting of the exponential form close to its
ends reveals the gap sizes, and also an indirect gap between these two momenta
forms at values of $t_{\perp}/t=4.5$, as reported in the main text. Compiling
these results, across a range of interplane hybridization values at $U/t=10$
[Fig. 6(c)], leads to a $\Delta_{\rm sp}$ largely consistent with the one
originally presented in Fig. 4(b), with a larger size dependence, however. The
minimum at $t_{\perp}/t\simeq 3$ still locates the transition from a
paramagnetic Mott insulator to the band-insulating regime.

##  Appendix D Analysis of a possible bond-order wave

The energy gaps displayed in Fig. 3(b) in the main text gave us indications
for the existence of multiple phases with growing $t_{\perp}$. In particular,
a possible candidate to explain the intervening regime in between the AFM Mott
insulator and the band-insulator at large interplane hybridization is a rung
bond-ordered wave (BOW). To account for this possibility, we calculate the BOW
order parameter,

$\langle\hat{B}_{i}\rangle\equiv\sum_{\sigma}(\hat{c}^{\dagger}_{i0\sigma}\hat{c}^{\phantom{\dagger}}_{i1\sigma}+\hat{c}^{\dagger}_{i1\sigma}\hat{c}^{\phantom{\dagger}}_{i0\sigma}),$
(4)

where, in the notation presented in the main text, the second sub-index refers
to the fermionic operator layer. Verification of the existence of long-range
order is accomplished via monitoring the BOW structure factor, $S_{\rm
BOW}=(1/L^{2})\sum_{ij}e^{-{\rm i}{\bf k}\cdot({\bf r}_{i}-{\bf
r}_{j})}\langle\hat{B}_{i}\hat{B}_{j}\rangle$.

In Fig. 7, the structure factors are shown as a function of the interlayer
hopping for two system sizes, as obtained for a system with $U/t=10$. In all
channels we investigate, ${\bf k}=[\pi,\pi],[0,\pi]$ or $[0,0]$, no apparent
size dependence of the BOW structure factor can be observed for a large range
of interplane hybridizations, suggesting the absence of such ordering in the
thermodynamic limit.

##  Appendix E Local correlators

While BOW rung correlations are short-ranged, other local correlators aid in
drawing a picture of the low-energy physics as the interplane hopping is
increased. The density-density correlations within a unit cell [Fig. 8(a)]
start at the value 0.25 when the planes are uncoupled ($t_{\perp}=0$) and
evolve such that equal-spin correlations are quickly suppressed whereas
opposite spin ones are enhanced until $t_{\perp}/t\simeq 2.1$. Accompanied by
the large negative interplane, intra-unit cell spin-spin correlations [Fig.
8(c)], these point to a robust interplane singlet formation in this regime, as
a local signature of the overall antiferromagnetic state. Past this threshold
in the hybridization $t_{\perp}$,
$\langle\hat{S}_{i,0}^{z}\hat{S}_{i,1}^{z}\rangle$ decreases in magnitude,
signaling the long-range order is absent and interplane singlets in different
unit cells become more and more independent. This can be seen explicitly in
Figs. 8(b) and 8(d), which report the interplane two-point correlations
(density and spin, respectively) for nearest-neighbor unit cells: they
asymptotically approach their uncorrelated values, pointing to the decreasing
interdependency of the different singlets within a unit cell. This approach is
continuous (at the low, but finite temperatures $T/t=1/20$) and does not
indicate a sharp transition to the collective product state of singlet states.
Rather, they are suggestive of a crossover, which is behind the explanation of
the smooth transition from the correlated to uncorrelated band-insulating
phases.

Figure 7: BOW structure factor $S_{\rm BOW}$ as a function of the interlayer
hopping $t_{\perp}/t$ for two system sizes at $U/t=10$ and $T/t=1/20$.
Different symmetry channels are displayed, (a) ${\bf k=[\pi,\pi]}$, (b)
$[\pi,0]$, and (c) $[0,0]$. The absence of any significative size dependence
indicates there is no long-range bond order.  Figure 8:  Top: Cartoon
depicting the different sites $(i,l)$ where the two-point correlations are
computed as displayed in the lower panels. Bottom: The interlayer (a),(b)
density-density correlations for same ($\sigma=\sigma^{\prime}$) and different
spin components ($\sigma\neq\sigma^{\prime}=-\sigma$) within the same unit-
cell (a) and nearest-neighbor unit cells (b). (c) and (d) display the
equivalent for the spin-spin correlations. All data are extracted at $U/t=10$
for $L=14$ and $T/t=1/20$. Spin summation is implicitly assumed for the
density correlations.

##  Appendix F Hamming Distance

Recent investigations in various fermionic models have shown that the onset of
criticality can be directly tracked by quantities related to metrics of the
auxiliary HS field [35, 50], bypassing the necessity of extraction of physical
observables. This field, local to each orbital in the spinful Hubbard model in
its $(d+1)$-dimensions formulation [22], is the quantity being sampled via a
Metropolis algorithm over the course of the Monte Carlo sampling [23]. A
useful metric (although others exist [51, 52]) is given by the average
distance traveled with respect to a given point in the phase space of
configurations $\\{s_{i,\tau}\\}$, and can be defined by the $L_{1}$-distance
(or Hamming distance) of configurations,

${\cal HD}=\frac{1}{2\cal V}\sum_{i,\tau}|s_{i,\tau}-s_{i,\tau}^{\rm ref.}|,$
(5)

with respect to a reference configuration $\\{s_{i,\tau}^{\rm ref.}\\}$. Here
${\cal V}$ is the length of the field, $L_{\tau}\cdot 2L^{2}$ in our case,
with $L_{\tau}$ being the number of imaginary-time slices that discretize the
inverse temperature: $\beta/L_{\tau}=\Delta\tau$. Assuming that
$\\{s_{i,\tau}^{\rm ref.}\\}$ is a typical configuration (namely, one obtained
after a significant number of warmup sweeps in the field), storing ${\cal HD}$
after each sweep on the sampling allows one to quantify the average distance
the Markov chain probes in the $(d+1)$-dimensional phase space. This was shown
to be intrinsically related to the physics of the models being investigated
[35]. For example, unordered phases were demonstrated to be associated with
completely uncorrelated HS configurations, rendering an average Hamming
distance, $\overline{\cal HD}$, equal to 0.5. Departure from this value
signals correlation within the sampled space, and thus physically ordered
phases.

Applying this idea to the bilayer Hubbard model, we report in Fig. 9 the
average Hamming distance $\overline{\cal HD}$ as a function of the interplane
hybridization at large interactions $U/t=10$. A significant departure from the
uncorrelated sampling occurs at $t_{\perp}\approx 5.8t$, which marks the
regime where the uncorrelated band-insulator gives way to a correlated one
[see Fig. 3(b) in the main text]. Another metric of the sampling is the
typical width $\sigma_{\cal HD}$ of the Gaussian-distributed Hamming distances
within a given realization. Averaging among independently seeded Markov
chains, $\overline{\sigma}_{\cal HD}$, results in a marked location for the
change of the distributions, which coincides with the onset of the ordered
phase at $t_{\perp}/t\approx 2.1$. That is, Hamming distances obtained at each
sampling process are much more diverse within the ordered phase in comparison
to the ones within the physically unordered regime. These results lend extra
insight about the location of the different phases, and have been checked to
exhibit qualitatively small finite-size effects.

Figure 9: (a) Average Hamming distance and averaged standard deviation of
Hamming distances versus interplane hybridization with $U/t=10$ and $L=14$.
(c) The corresponding spin-resolved sign $\langle{\cal S}_{\sigma}\rangle$.
Vertical shaded regions in (a) and (d) display the QCP obtained by AFM
structure factor $t^{c}_{\perp}/t\approx 2.1$, whereas the shaded regions in
the zoomed-in panels (b) and (d) mark the confidence region in $t_{\perp}$
that indicates the correlated-to-uncorrelated band insulating crossover (see
main text).

##  Appendix G ED results

We start by analyzing the strongly interacting regimes, where finite-size
effects are typically less dramatic: Figs. 10(a)-10(d) for $U/t=10$, and
Figs.10(e)-10(h) for $U/t=32$. The dependence of the antiferromagnetic spin
structure factor (see main text for definition) is shown in Fig.10(a)
and10(e). As the system evolves from a typical planar antiferromagnet to a
bilayer one as $t_{\perp}$ increases, i.e., the number of neighbors
effectively grows, this quantity (computed in the ground-state
$|\Psi_{0}\rangle$) initially increases. However, at larger interplane
hybridization, it is then suppressed for $t_{\perp}\gtrsim t$ as singlets
form, similar to the behavior seen in the DQMC results in the main text.
Concomitantly, given the reduced spin ordering, the staggered charge structure
factor, $S_{\rm
cdw}\equiv(1/2L^{2})\sum_{i,j}(-1)^{i+j}\langle(\hat{n}_{i\uparrow}+\hat{n}_{i\downarrow})(\hat{n}_{j\uparrow}+\hat{n}_{j\downarrow})\rangle$
slightly increases with $t_{\perp}$, albeit with an overall small magnitude,
suggesting the absence of a charge density wave formation. The smooth
evolution of these two quantities is evidence against the manifestation of a
first-order phase transition in the range of parameters studied.

Figure 10: ED results of observables in a $\sqrt{8}\times\sqrt{8}$ bilayer at
$U/t=10$ versus interplane hybridization: (a) the expectation value of the
staggered charge and spin-structure factor, (b) intra- and interplane spin-
spin nearest neighbor correlations, (c) double occupancy, and (d) the fidelity
susceptibility. (e)-(h) The same for $U=32t$; vertical dashed line in (h)
depicts the Heisenberg result from Refs. [2, 20] $J^{c}_{\perp}/J=2.52181(3)$
converted to $t^{c}_{\perp}/t=1.58802(1)$ via their correspondence in the
strongly interacting regime, $J\propto t^{2}/U$.

The reduction of global antiferromagnetic order giving way to a unordered
state displaying interplane spin singlets can be also seen by directly
computing the nearest-neighbor intra- and interplane spin correlation
functions,
$\langle\hat{m}_{i}\hat{m}_{j}\rangle\equiv\langle(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow})(\hat{n}_{j\uparrow}-\hat{n}_{j\downarrow})\rangle$.
All such spin correlations are negative, which signals their local
antiferromagnetic character. With increasing $t_{\perp}$, interplane spin
correlations surpass intraplane ones in magnitude, preceding the onset of the
quantum disordered phase.

The double occupancy,
$\langle\hat{n}_{\uparrow\downarrow}\rangle\equiv\langle\hat{n}_{\uparrow}\hat{n}_{\downarrow}\rangle$,
in analogy with Figs. 4(b) and 4(d), can give the locations of the insulator-
metallic pseudo-(or finite-size influenced) transitions at small interactions,
in particular, when its derivative,
$d\langle\hat{n}_{\uparrow\downarrow}\rangle/dt_{\perp}$, is evaluated. At
large $U/t$, Fig. 10(g), their absolute value is relatively small, and
reflects the approach to the Heisenberg regime of localized spins.

While much can be captured by the dependence of few-body correlators as we
have done so far, a direct account of how the many-body ground state changes
with the increasing interplane hybridization, can be understood in terms of
the fidelity susceptibility
$g_{\perp}=\frac{1}{L^{2}}\frac{1-|\langle\Psi_{0}(t_{\perp})|\Psi_{0}(t_{\perp}+dt_{\perp})\rangle|}{dt_{\perp}^{2}}$
[53, 54, 55, 56], which can locate a quantum phase transition without making
direct assumptions regarding a possible order parameter [57, 58, 59]. This
quantity displays a peak that is extensive, and independent of the
(sufficiently small) parameter ‘perturbation’ for continuous phase
transitions. Here, we take $dt_{\perp}=10^{-3}t$, and we note the presence of
a first peak at small $t_{\perp}$ as related to the crossover from planar to
bilayer antiferromagnetism. A second peak is also visible, which in the case
of almost fully formed local spins captures the magnetic transition, closely
following the results obtained for the Heisenberg model for much larger
lattices [2, 20].

A more general description of such quantities across the $t_{\perp}-U$ phase
diagram is given in Fig. 11. At small interactions, the various observables
are deeply affected by finite-size effects, manifesting the influence of the
metallic region extending at finite values of $U$. In this small cluster, only
eight $k$ points are available in the Brillouin zone. In this case, the
nesting condition (see main text) between the two bands occurs only at
$t_{\perp}/t=0$ and 4. As a result, the phase diagrams are remarkably similar
to the ones obtained in small clusters in DMFT [17]. In particular, we can
highlight (i) the large antiferromagnetic structure factor at small
$t_{\perp}$ denoting the magnetically ordered phase [Fig. 11(a)]; (ii) the
$d\langle\hat{n}_{\uparrow\downarrow}\rangle/dt_{\perp}>0$
[$d\langle\hat{n}_{\uparrow\downarrow}\rangle/dt_{\perp}<0$] describing the
insulator-metal [metal-insulator] transition for this small cluster [Fig.
11(b)], similar to the results of Fig. 4; (iii) the onset of the quantum
disordered phase, signified by the close-to-saturation of the interplane
nearest-neighbor spin correlations, $\langle\hat{m}_{i}\hat{m}_{j}\rangle_{\rm
inter}$ in Fig. 11(c); and (iv) the second peak (branch) of the fidelity
susceptibility at $t_{\perp}/t\approx 4$ describes the magnetic transition at
small interactions, and converges to the Heisenberg limit at sufficiently
large values of $U/t$ [Fig. 11(d)].

Figure 11: Color maps of several quantities in the space of parameters
$t_{\perp}$ versus $U$ for the $2\times 8$-site bilayer lattice. (a)
Expectation value of the antiferromagnetic spin-structure factor, (b)
derivative of the double occupancy in respect to $t_{\perp}$, (c) interplane
spin-spin nearest-neighbor correlations, and (d) the fidelity susceptibility
$g_{t_{\perp}}$. Star markers depict the magnetic transition at $T=0$ obtained
using QMC in much larger lattices in Ref. [15]; cross marker on the right
(large $U$) panel edge denotes the Heisenberg limit [2, 20]; and diamond-
shaped ones mark the nesting conditions in the noninteracting regime for this
cluster size.

Lastly, we provide a simplified view of the single-particle gap within ED
calculations by computing the charge-gap for excitations,

$\Delta_{c}=[E_{0}(N_{e}+1)-E_{0}(N_{e})]-[E_{0}(N_{e})-E_{0}(N_{e}-1)],$ (6)

where $E_{0}(N_{e})$ is the ground-state energy with $N_{e}$ electrons. Since
the lattice is small, this quantity displays relatively large finite-size
effects, but in principle should capture similar information as the single-
particle gap extracted from QMC simulations in the main text, in particular,
in the case where twisted-boundary conditions are applied [60].

Here, using standard periodic boundary conditions, we report $\Delta_{c}$ in
Fig. 12, accompanied by the corresponding fidelity susceptibility. We avoid
the regime of small interactions, focusing on $U/t=10$, and contrast it with
the results in much larger systems obtained via QMC in Fig. 3(b) of the main
text.

Figure 12: The charge gap $\Delta_{c}$ (right-axis) versus interplane
hybridization for a $\sqrt{8}\times\sqrt{8}$ bilayer, computed via ED with
$U/t=10$. For comparison, we also show the fidelity susceptibility (left-axis)
at the same interaction strength. The local minimum (maximum) of $\Delta_{c}$
($g_{t_{\perp}}$) is marked via a dotted (dashed) vertical line at
$t_{\perp}/t=3.3$ ($t_{\perp}/t=2.8$).

Although at small values of $t_{\perp}$ the gaps $\Delta_{c}$ and $\Delta_{\rm
sp}$, extracted from ED and QMC respectively, do not display similar behavior
(notwithstanding the strikingly different lattice sizes), at large
hybridization, both exhibit a minimum at around $t_{\perp}/t\simeq 3$, which
we have associated to a crossover from a paramagnetic Mott insulator towards a
band insulator in the main text. To argue that the magnetic transition is not
aligned with the above described crossover, we overlay the data originally
shown in Fig. 10(d) for the fidelity susceptibility in Fig. 12. The peak at
large $t_{\perp}/t$, which we inferred to be related to the magnetic
transition, does not coincide with the local minimum of $\Delta_{c}$, opening
room for a potential paramagnetic Mott insulator surviving at
$2.8<t_{\perp}/t<3.3$ for this lattice size and for $U/t=10$.

## References

  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Possible high Tc superconductivity in the Ba-La-Cu-O system, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Sandvik and Scalapino [1994] A. W. Sandvik and D. J. Scalapino, Order-disorder transition in a two-layer quantum antiferromagnet, Phys. Rev. Lett. 72, 2777 (1994).
  * Sandvik and Scalapino [1996] A. W. Sandvik and D. J. Scalapino, Effects of intrabilayer coupling on the magnetic properties of YBa2Cu3O6, Phys. Rev. B 53, R526 (1996).
  * Lercher and Wheatley [1994] M. J. Lercher and J. M. Wheatley, Magnetic bilayer coupling in YBa2Cu3O7, Phys. Rev. B 49, 736 (1994).
  * Eder  _et al._ [1995] R. Eder, Y. Ohta, and S. Maekawa, Ground-state properties and dynamics of the bilayer $t$-${J}$ model, Phys. Rev. B 52, 7708 (1995).
  * Bulut  _et al._ [1992] N. Bulut, D. J. Scalapino, and R. T. Scalettar, Nodeless d-wave pairing in a two-layer Hubbard model, Phys. Rev. B 45, 5577 (1992).
  * Scalettar  _et al._ [1994] R. T. Scalettar, J. W. Cannon, D. J. Scalapino, and R. L. Sugar, Magnetic and pairing correlations in coupled Hubbard planes, Phys. Rev. B 50, 13419 (1994).
  * dos Santos [1995] R. R. dos Santos, Magnetism and pairing in Hubbard bilayers, Phys. Rev. B 51, 15540 (1995).
  * Loh  _et al._ [1990] E. Loh, J. Gubernatis, R. Scalettar, S. White, D. Scalapino, and R. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Quantum critical points and the sign problem, Science 375, 418 (2022).
  * Rademaker  _et al._ [2013] L. Rademaker, S. Johnston, J. Zaanen, and J. van den Brink, Determinant quantum Monte Carlo study of exciton condensation in the bilayer Hubbard model, Phys. Rev. B 88, 235115 (2013).
  * Bouadim  _et al._ [2008] K. Bouadim, G. G. Batrouni, F. Hébert, and R. T. Scalettar, Magnetic and transport properties of a coupled Hubbard bilayer with electron and hole doping, Phys. Rev. B 77, 144527 (2008).
  * Euverte  _et al._ [2013] A. Euverte, S. Chiesa, R. T. Scalettar, and G. G. Batrouni, Magnetic transition in a correlated band insulator, Phys. Rev. B 87, 125141 (2013).
  * Golor  _et al._ [2014] M. Golor, T. Reckling, L. Classen, M. M. Scherer, and S. Wessel, Ground-state phase diagram of the half-filled bilayer Hubbard model, Phys. Rev. B 90, 195131 (2014).
  * Fuhrmann  _et al._ [2006] A. Fuhrmann, D. Heilmann, and H. Monien, From Mott insulator to band insulator: A dynamical mean-field theory study, Phys. Rev. B 73, 245118 (2006).
  * Kancharla and Okamoto [2007] S. S. Kancharla and S. Okamoto, Band insulator to Mott insulator transition in a bilayer Hubbard model, Phys. Rev. B 75, 193103 (2007).
  * Hafermann  _et al._ [2009] H. Hafermann, M. I. Katsnelson, and A. I. Lichtenstein, Metal-insulator transition by suppression of spin fluctuations, EPL (Europhysics Letters) 85, 37006 (2009).
  * Ruger  _et al._ [2014] R. Ruger, L. F. Tocchio, R. Valenti, and C. Gros, The phase diagram of the square lattice bilayer Hubbard model: A variational Monte Carlo study, New Journal of Physics 16, 033010 (2014).
  * Wang  _et al._ [2006] L. Wang, K. S. D. Beach, and A. W. Sandvik, High-precision finite-size scaling analysis of the quantum-critical point of $s=1/2$ Heisenberg antiferromagnetic bilayers, Phys. Rev. B 73, 014431 (2006).
  * Karakuzu  _et al._ [2021] S. Karakuzu, S. Johnston, and T. A. Maier, Superconductivity in the bilayer Hubbard model: Are two Fermi surfaces better than one? (2021), arXiv:2107.13996 .
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte Carlo calculations of coupled boson-fermion systems. I, Phys. Rev. D 24, 2278 (1981).
  * Hirsch [1985] J. E. Hirsch, Two-dimensional Hubbard model: Numerical simulation study, Phys. Rev. B 31, 4403 (1985).
  * [24] R. Mondaini, S. Tarat, and R. Scalettar, submitted.
  * White  _et al._ [1989] S. White, D. Scalapino, R. Sugar, E. Loh, J. Gubernatis, and R. Scalettar, Numerical study of the two-dimensional Hubbard model, Phys. Rev. B 40, 506 (1989).
  * Paiva  _et al._ [2005] T. Paiva, R. Scalettar, W. Zheng, R. Singh, and J. Oitmaa, Ground-state and finite-temperature signatures of quantum phase transitions in the half-filled Hubbard model on a honeycomb lattice, Phys. Rev. B 72, 085123 (2005).
  * Meng  _et al._ [2010] Z. Meng, S. Wessel, A. Muramatsu, T. Lang, and F. Assaad, Quantum spin liquid emerging in two-dimensional correlated Dirac fermions, Nature 464, 847 (2010).
  * Sorella  _et al._ [2012] S. Sorella, Y. Otsuka, and S. Yunoki, Absence of a spin liquid phase in the Hubbard model on the honeycomb lattice, Sci. Rep. 2, 992 (2012).
  * Note [1] Exact diagonalization results in Appendix G provide further evidence that the magnetic and the Mott insulating transitions are not aligned at large interaction strengths in the bilayer Hubbard model.
  * Wang  _et al._ [2020] Y. D. Wang, W. L. Yao, Z. M. Xin, T. T. Han, Z. G. Wang, L. Chen, C. Cai, Y. Li, and Y. Zhang, Band insulator to Mott insulator transition in 1T-TaS2, Nature Communications 11, 4215 (2020).
  * Sandvik [1998] A. W. Sandvik, Stochastic method for analytic continuation of quantum Monte Carlo data, Phys. Rev. B 57, 10287 (1998).
  * Note [2] We are currently unaware of other possible orders that may manifest in the described paramagnetic Mott insulator. A rung bond-ordered wave is not a good candidate – see Appendix D.
  * Lee  _et al._ [2014] H. Lee, Y.-Z. Zhang, H. O. Jeschke, and R. Valentí, Competition between band and Mott insulators in the bilayer Hubbard model: A dynamical cluster approximation study, Phys. Rev. B 89, 035139 (2014).
  * Sentef  _et al._ [2009] M. Sentef, J. Kuneš, P. Werner, and A. P. Kampf, Correlations in a band insulator, Phys. Rev. B 80, 155116 (2009).
  * Yi  _et al._ [2021] T.-C. Yi, R. T. Scalettar, and R. Mondaini, Hamming distance and the onset of quantum criticality (2021).
  * Garg  _et al._ [2006] A. Garg, H. Krishnamurthy, and M. Randeria, Can correlations drive a band insulator metallic?, Phys. Rev. Lett. 97, 046403 (2006).
  * Paris  _et al._ [2007] N. Paris, K. Bouadim, F. Hebert, G. Batrouni, and R. Scalettar, Quantum Monte Carlo study of an interaction-driven band-insulator–to–metal transition, Phys. Rev. Lett. 98, 046403 (2007).
  * Bouadim  _et al._ [2007] K. Bouadim, N. Paris, F. Hébert, G. Batrouni, and R. Scalettar, Metallic phase in the two-dimensional ionic Hubbard model, Phys. Rev. B 76, 085112 (2007).
  * Craco  _et al._ [2008] L. Craco, P. Lombardo, R. Hayn, G. Japaridze, and E. Müller-Hartmann, Electronic phase transitions in the half-filled ionic Hubbard model, Phys. Rev. B 78, 075121 (2008).
  * Garg  _et al._ [2014] A. Garg, H. Krishnamurthy, and M. Randeria, Doping a correlated band insulator: A wew route to half-metallic behavior, Phys. Rev. Lett. 112, 106406 (2014).
  * Koepsell  _et al._ [2020] J. Koepsell, S. Hirthe, D. Bourgund, P. Sompet, J. Vijayan, G. Salomon, C. Gross, and I. Bloch, Robust bilayer charge pumping for spin- and density-resolved quantum gas microscopy, Phys. Rev. Lett. 125, 010403 (2020).
  * Gall  _et al._ [2021] M. Gall, N. Wurz, J. Samland, C. F. Chan, and M. Köhl, Competing magnetic orders in a bilayer Hubbard model with ultracold atoms, Nature 589, 40 (2021).
  * Feng  _et al._ [2001] D. L. Feng, N. P. Armitage, D. H. Lu, A. Damascelli, J. P. Hu, P. Bogdanov, A. Lanzara, F. Ronning, K. M. Shen, H. Eisaki, C. Kim, Z.-X. Shen, J.-i. Shimoyama, and K. Kishio, Bilayer splitting in the electronic structure of heavily overdoped ${Bi}_{2}{Sr}_{2}{CaCu}_{2}{O}_{8+\delta}$, Phys. Rev. Lett. 86, 5550 (2001).
  * Liechtenstein  _et al._ [1996] A. I. Liechtenstein, O. Gunnarsson, O. K. Andersen, and R. M. Martin, Quasiparticle bands and superconductivity in bilayer cuprates, Phys. Rev. B 54, 12505 (1996).
  * Lanatà  _et al._ [2009] N. Lanatà, P. Barone, and M. Fabrizio, Superconductivity in the doped bilayer Hubbard model, Phys. Rev. B 80, 224524 (2009).
  * Maier and Scalapino [2011] T. A. Maier and D. J. Scalapino, Pair structure and the pairing interaction in a bilayer Hubbard model for unconventional superconductivity, Phys. Rev. B 84, 180513 (2011).
  * Matsumoto  _et al._ [2020] K. Matsumoto, D. Ogura, and K. Kuroki, Strongly enhanced superconductivity due to finite energy spin fluctuations induced by an incipient band: A FLEX study on the bilayer Hubbard model with vertical and diagonal interlayer hoppings, Journal of the Physical Society of Japan 89, 044709 (2020).
  * Bohrdt  _et al._ [2022] A. Bohrdt, L. Homeier, I. Bloch, E. Demler, and F. Grusdt, Strong pairing in mixed-dimensional bilayer antiferromagnetic Mott insulators, Nature Physics 10.1038/s41567-022-01561-8 (2022).
  * [49] QUEST: QUantum Electron Simulation Toolbox, http://www.cs.ucdavis.edu/bai/QUEST/.
  * Tirelli and Costa [2021] A. Tirelli and N. C. Costa, Learning quantum phase transitions through topological data analysis, Phys. Rev. B 104, 235146 (2021).
  * Mendes-Santos  _et al._ [2021a] T. Mendes-Santos, X. Turkeshi, M. Dalmonte, and A. Rodriguez, Unsupervised learning universal critical behavior via the intrinsic dimension, Phys. Rev. X 11, 011040 (2021a).
  * Mendes-Santos  _et al._ [2021b] T. Mendes-Santos, A. Angelone, A. Rodriguez, R. Fazio, and M. Dalmonte, Intrinsic dimension of path integrals: Data-mining quantum criticality and emergent simplicity, PRX Quantum 2, 030332 (2021b).
  * Zanardi and Paunković [2006] P. Zanardi and N. Paunković, Ground state overlap and quantum phase transitions, Phys. Rev. E 74, 031123 (2006).
  * Campos Venuti and Zanardi [2007] L. Campos Venuti and P. Zanardi, Quantum critical scaling of the geometric tensors, Phys. Rev. Lett. 99, 095701 (2007).
  * Zanardi  _et al._ [2007] P. Zanardi, P. Giorda, and M. Cozzini, Information-theoretic differential geometry of quantum phase transitions, Phys. Rev. Lett. 99, 100603 (2007).
  * You  _et al._ [2007] W.-L. You, Y.-W. Li, and S.-J. Gu, Fidelity, dynamic structure factor, and susceptibility in critical phenomena, Phys. Rev. E 76, 022101 (2007).
  * Yang [2007] M.-F. Yang, Ground-state fidelity in one-dimensional gapless models, Phys. Rev. B 76, 180403 (2007).
  * Jia  _et al._ [2011] C. J. Jia, B. Moritz, C.-C. Chen, B. S. Shastry, and T. P. Devereaux, Fidelity study of the superconducting phase diagram in the two-dimensional single-band Hubbard model, Phys. Rev. B 84, 125113 (2011).
  * Mondaini  _et al._ [2015] R. Mondaini, P. Nikolić, and M. Rigol, Mott-insulator–to–superconductor transition in a two-dimensional superlattice, Phys. Rev. A 92, 013601 (2015).
  * Poilblanc [1991] D. Poilblanc, Twisted boundary conditions in cluster calculations of the optical conductivity in two-dimensional lattice models, Phys. Rev. B 44, 9562 (1991).


0 Research  Our visualization system with three coordinated views. The Tree
view (a) presents the structure of an embedding tree and data distributions
across branches. The Dimensionality Reduction view (b) shows data entities
from a selected tree node and reveals the cluster information. The Data Table
(c) shows details of the selected data entities and enables users to
search/filter/sort them.

#  EmbeddingTree: Hierarchical Exploration of Entity Features in Embedding

Yan Zheng e-mail<EMAIL_ADDRESS>Junpeng Wang<EMAIL_ADDRESS>Chin-
Chia Michael Yeh<EMAIL_ADDRESS>Yujie Fan<EMAIL_ADDRESS>Huiyuan
Chen<EMAIL_ADDRESS>Liang Wang<EMAIL_ADDRESS>Wei Zhang
<EMAIL_ADDRESS>Visa Research

###### Abstract

Embedding learning transforms discrete data entities into continuous numerical
representations, encoding features/properties of the entities. Despite the
outstanding performance reported from different embedding learning algorithms,
few efforts were devoted to structurally interpreting how features are encoded
in the learned embedding space. This work proposes EmbeddingTree, a
hierarchical embedding exploration algorithm that relates the semantics of
entity features with the less-interpretable embedding vectors. An interactive
visualization tool is also developed based on EmbeddingTree to explore high-
dimensional embeddings. The tool helps users discover nuance features of data
entities, perform feature denoising/injecting in embedding training, and
generate embeddings for unseen entities. We demonstrate the efficacy of
EmbeddingTree and our visualization tool through embeddings generated for
industry-scale merchant data and the public 30Music listening/playlists
dataset.

##  1 Introduction

Embedding learning has emerged as a powerful technique for analyzing text
documents [1, 2] and other data types [3, 4, 5, 6, 7]. It encapsulates
relevant information into high-dimensional (HD) vectors that could be utilized
for various downstream tasks. Often, there is no clear correspondence between
raw features and individual embedding dimensions, and each dimension could
entangle the semantics of multiple features. This flexibility allows for
embedding more features and better capturing the correlations between them.
However, it also eschews the interpretability of the embeddings.

Targeting on the interpretability problem of embeddings, multiple visual
exploration solutions have been proposed [8, 9, 10, 11, 12]. These solutions
explore embeddings by relating their coordinates’ similarity to specific
features or feature combinations. For example, a cluster of embeddings with
similar spatial locations in embedding projector [8] often corresponds to data
entities with similar feature values. The correspondence helps to disclose the
encoded semantics in the embedding. We consider this type of exploration
parallel exploration, as features are treated equally and explored in
parallel.

In some datasets, however, features are not of the same importance and should
be explored in order or following a hierarchy. For example, in the merchant
embeddings generated by Wang et al. [13] using credit-card transactions,
merchants’ location is the feature that dominates the embedding the most,
visiting frequency is the second, and merchants’ category code (MCC) is the
third. As shown in Figure 1a, when using PCA to project the merchant
embeddings, the merchants get separated into four clear clusters, representing
merchants from four different cities. Focusing on the red and cyan clusters
(merchants from Los Angeles and Manhattan), MCC and visiting frequency
dominate the layouts of the two local clusters respectively (Figure 1b, c).
Existing parallel explorations cannot easily reveal this hierarchical feature
dominance and could hinder the exploration process from discovering feature-
level insights.

To accommodate the aforementioned problem, we introduce hierarchical
exploration in this work, where features form nested structures, and users can
explore them in a top-down hierarchy. The hierarchy in Figure 1 could be
extracted as a tree structure, as highlighted in the light yellow region of
Figure EmbeddingTree: Hierarchical Exploration of Entity Features in
Embeddinga, guiding the exploration process. The benefit of such exploration
is twofold: (1) the hierarchy directly infers features’ importance in the
embedding; (2) the hierarchy provides explicit guidance and helps users drill
down to embeddings with interested features more efficiently.

Figure 1: Hierarchical structure in the embedding. (a) merchant embeddings are
dominated by the location feature and separated into four clusters for the
four cities. (b, c) For the merchants from Los Angeles (LA) and Manhattan
(MAN), MCC and visiting frequency further dominate the embedding layouts,
respectively. The tree branches highlighted in light-yellow in Figure
EmbeddingTree: Hierarchical Exploration of Entity Features in Embeddinga show
this hierarchy.

In detail, we first propose a Gaussian Mixture Model (GMM) based algorithm to
generate an EmbeddingTree from HD embeddings. The algorithm uses entity
features to split the embedding space based on which feature can best separate
the embedding in nested local regions through Bayesian Information Criterion
(BIC). Next, we develop a visualization tool (Figure EmbeddingTree:
Hierarchical Exploration of Entity Features in Embedding) to visualize the
trained embeddings and organize them according to the hierarchical entity
features in the EmbeddingTree. The tool employs interactive tree
visualization, dimensionality reduction, and data table visualization to
facilitate the hierarchical exploration of HD embeddings. Our EmbeddingTree
algorithm and the visual exploration can help to analyze the hierarchical
importance of entity features in the embeddings. As a result, users could
discover the inconsistency between entity features and embeddings, as well as
generate the embeddings for new entities with little structure information but
rich entity features.

To sum up, our contributions in this work are:

  1. 1.

We propose a GMM-based algorithm to extract a feature hierarchy from HD
embeddings.

  2. 2.

We introduce a visual analytics tool that helps users effectively explore
embedding data following the extracted hierarchy.

##  2 Related Work

#### Embedding Generation

The embedding technology is popularized by the introduction of Word2Vec [14]
and GloVe [15]. Both were proposed for words and rely on the distributional
hypothesis of Harris [16], i.e., words occurring in similar contexts tend to
have similar meanings. Word embedding generated from those methods captures
the semantic and syntactic properties of words. Another popular group of
embeddings are graph embeddings, which aim to represent graphs or graph nodes
as HD vectors/embeddings. These embeddings capture the essential information
from both the structure-level and node-level of graphs and help to boost the
downstream tasks’ performance, e.g., node classification [17], recommendation
[18] and link prediction [19]. With the fast evolution of deep learning
techniques, more and more embedding learning techniques [20, 21, 22, 23, 24]
have been proposed, and they have been used to generate embeddings for
different types of data.

Our EmbedingTree algorithm and the visualization tool can be applied to the
embedding dataset with categorical/numerical features, regardless of the
algorithm to generate the embeddings. For example, the merchant embeddings in
Section 5.1 were trained from credit-card transaction sequences, using the the
Word2vec [25] algorithm; the user embeddings in Section 5.2 were trained from
the 30Music dataset, using the embedding learning framework proposed in [26].

#### Embedding Visualization

Existing embedding visualization works can roughly be categorized into three
groups, focusing on semantic interpretation, analogy reasoning, and comparison
of embeddings. For the first group, Latent Space Cartography [10] and
LatentVis [27] disclose how semantics are encoded and expressed along certain
directions of the HD latent/embedding space. For the second group, Liu et al.
[11] study the analogy relation between pairs of word embeddings and optimize
embedding generation algorithms by better preserving the relation. For the
third group of embedding comparisons, EmbeddingVis [28] compares graph
embeddings generated from different methods by linking multiple juxtaposed
scatterplots. embComp [9] quantifies the similarity of two embeddings through
the number of shared local neighbors. Recently, Embedding Comparator[29]
presents a global comparison of embedding spaces alongside fine-grained
inspection of local neighborhoods.

Most of these embedding visualizations explore embedding data without
explicitly leveraging the relationship between entity features. Our work
targets at this missing piece and strives to explore embeddings by following a
hierarchy of data features distilled from the embedding data. This helps to
structurally explore the embedding space and better understand the encoded
semantics.

##  3 Algorithm

Our EmbeddingTree construction is based on the traditional decision tree
algorithm but uses two sets of data, one is the embedding data and the other
is entity features. Entity features are used as the search criteria to split
the embedding by efficiently computing an approximate Bayesian Information
Criterion (BIC) for a Gaussian Mixture Model (GMM).

###  3.1 Embedding Tree Algorithm

The input of EmbeddingTree is a combination of an embedding dataset
$\mathbf{X}=\\{\mathbf{x}_{1},...,\mathbf{x}_{N}\\}$
($\mathbf{x}_{i}\in\mathbb{R}^{p}$, $1\leq i\leq N$) and the corresponding
feature set. For numerical and multi-categorical feature, we pre-process the
features by (1) binning numerical features into categorical ones; (2)
converting non-binary categorical features into binary ones by asking yes-no
questions regarding to each of their categorical values. Finaly we have the
binary feature sets: $\mathbf{F}=\\{\mathbf{f}_{1},...,\mathbf{f}_{N}\\}$
($\mathbf{f}_{i}\in\\{0,1\\}^{q}$, $1\leq i\leq N$) where $N$ is the number of
entities and $\mathbf{f}_{i}$ is the feature set of the entity with embedding
$\mathbf{x}_{i}$. $p$ and $q$ denote the embeddings’ dimensionality and the
total number of converted binary features.

With the processed data, we describe the EmbeddingTree algorithm with details
in Algorithm 1. First, before a tree splitting, PCA is applied to the
embedding data to find the first eigenvector and project all embeddings onto
the principal direction. Next, we iterate through the $q$ features (line 6)
and evaluate them based on the splitting criteria described in Section A.1 to
pick out the best feature for splitting (line 8-10), using the feature’s
binary value (line 11-13). Lastly, the above procedure is executed recursively
(line 15-16) until the splitting criterion ($\Theta$), e.g., the number of
entities per tree node or the tree depth, is no longer satisfied (line 2).

With the given embedding and feature data, the whole procedure is
deterministic. The goal of constructing EmbeddingTree is to provide users a
tool to hierarchically explore and understand the embeddings by incorporating
the entity features.

Algorithm 1 Build an EmbeddingTree

1:procedure BuildTree($[\mathbf{X},\mathbf{F}],q,\Theta$)

2: if $\Theta\hskip 2.84526ptis\hskip 2.84526ptnot\hskip 2.84526ptsatisfied$
then

3: return LeafNode([$\mathbf{X},\mathbf{F}$])

4: else

5: $max\\_t\leftarrow-\infty$

6: for $k\in\\{1,...,q\\}$ do

7: $t=Embedding\textendash BIC([\mathbf{X},\mathbf{F}^{k}])$

8: if $t>max\\_t$ then

9: $bestFea=k$

10: $max\\_t=t$

11:
$[\mathbf{X},\mathbf{F}]_{left}=\\{\mathbf{x}\in\mathbf{X}|\mathbf{F}_{bestFea}==0\\}$

12:

13:
$[\mathbf{X},\mathbf{F}]_{right}=\\{\mathbf{x}\in\mathbf{X}|\mathbf{F}_{bestFea}==1\\}$

14:

15: Children.Left = $BuildTree([\mathbf{X},\mathbf{F}]_{left},q,\Theta)$

16: Children.Right = $BuildTree([\mathbf{X},\mathbf{F}]_{right},q,\Theta)$

17: return Children

###  3.2 2-GMM Splitting with Embedding-BIC.

One critical component of Algorithm 1 is the criterion used to select the best
splitting feature. The criterion is computed based on the approximate BIC for
GMMs inspired by [30]. Specifically, each binary feature value splits the
embeddings into two clusters, each is then formulated as a Gaussian. Suppose
one feature $\mathbf{F}^{k}$, $k\in[1,q]$ splits the first $s$ entities with
feature value $0$ and the rest $N-s$ entities with feature value 1, the
complete likelihood is:

$P(x,\mu,\sigma,w)=\prod_{i=1}^{s}w_{1}\mathcal{N}(x_{i};\mu_{1},{\sigma_{1}}^{2})\prod_{i=s+1}^{N}w_{2}\mathcal{N}(x_{i};\mu_{2},{\sigma_{2}}^{2}),$
(1)

where $\mathcal{N}$ is the probability density function of Gaussian
distribution. For $i\in(1,2)$, $\mu_{i}$ is the mean, $\sigma_{i}$ is the
standard deviation of the Gaussian functions, $w_{i}$ is the weight of the
cluster and $x_{i}$ is the PCA projected value of the embedding
$\mathbf{x}_{i}$.

We estimate the means ($\hat{\mu}_{i}$), variances ($\hat{\sigma}_{i}$) and
weights ($\hat{w}_{i}$) for both clusters using the maximum likelihood
estimation (MLE).

$\hat{\mu_{1}}=\frac{1}{s}\sum_{i=1}^{s}x_{i},\hskip
14.22636pt\hat{\sigma_{1}}=\frac{1}{s}\sum_{i=1}^{s}||x_{i}-\hat{\mu_{1}}||^{2},\hskip
14.22636pt\hat{w_{1}}=\frac{s}{N},$
$\hat{\mu_{2}}=\frac{1}{N-s}\sum_{i=s+1}^{N}x_{i},\hskip
14.22636pt\hat{\sigma_{2}}=\frac{1}{N-s}\sum_{i=s+1}^{N}||x_{i}-\hat{\mu_{2}}||^{2},\hskip
14.22636pt\hat{w_{2}}=\frac{N-s}{N}.$

Plugging in the MLE for all the parameters, substituting the probability
density function of Gaussian
$\mathcal{N}=\frac{1}{\sqrt{2\pi}\sigma}\exp^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}$
into Equation 1 and simplifying the outcome, we get the following expression
for the log likelihood for any given $s$:

$\hat{L}_{s}=-\frac{s}{2}\log 2\pi\hat{\sigma_{1}}^{2}-\frac{N-s}{2}\log
2\pi\hat{\sigma_{2}}^{2}\\\ +s\log\hat{w_{1}}+(N-s)\log\hat{w_{2}}.$ (2)

We call our algorithm Embedding-BIC algorithm. Specifically, it chooses the
feature as the split-point that maximizes $\hat{L_{s}}$ . The detailed
explanation and pseudocode for this algorithm can be found in our Appendix.

###  3.3 Embedding and Entity Feature Consistency

EmbeddingTree splits the whole dataset (embedding data + feature data) into
different clusters represented by the leaf nodes. The relationship between
entity features and embeddings can be classified into two categories,
consistent and inconsistent.

If entity features and embeddings are consistent, the embeddings fully encode
the feature information, and they can be perfectly separated by the
EmbeddingTree, i.e., each leaf node contains a single cluster of embeddings.
In the second case, either the embedding doesn’t show the cluster information
under a subset of features or there are multiple clusters in one leaf node.
The user can dive into the problematic clusters to diagnose the reasons for
the inconsistency.

##  4 Visualization Tool

We develop a visualization system to help users interactively explore the
generated EmbeddingTree. The system contains three main visualization
components, as shown in Figure EmbeddingTree: Hierarchical Exploration of
Entity Features in Embedding.

The Tree view (Figure EmbeddingTree: Hierarchical Exploration of Entity
Features in Embeddinga) presents the structure of the EmbeddingTree. The texts
of internal nodes denote the splitting condition for each branch of the tree,
whereas the numbers in tree leaves represent the number of entities falling
into the corresponding leaves. The links connecting a parent and a child node
show the flowing of data entities (branched by different splitting
conditions), and their width is proportional to the number of entities.

For the exploration of large and deep trees, we empower the tree visualization
with interactive collapsing of the tree nodes. As shown in Figure
EmbeddingTree: Hierarchical Exploration of Entity Features in Embedding-a1,
clicking individual internal nodes will expand/fold the corresponding tree
branches to save the screen space (e.g., the light blue nodes have been
folded). Clicking leaf nodes will select the entities of the corresponding
tree branch into the Dimensionality Reduction view for detailed exploration.

The Dimensionality Reduction view (Figure EmbeddingTree: Hierarchical
Exploration of Entity Features in Embeddingb) demonstrates the selected data
entities by projecting their HD embedding into 2D for visualization and
interaction. The selected leaf in Figure EmbeddingTree: Hierarchical
Exploration of Entity Features in Embeddinga (highlighted with the red dashed
circle) has 471 entities, and they are represented by the 471 points in the
Dimensionality Reduction view. We have explored the three most popular
dimensionality reduction algorithms, i.e., tSNE, PCA, and UMAP. All meet our
requirement of disclosing the cluster pattern from the selected data entities.
We choose PCA as the default algorithm considering its significantly lower
computational cost. With the visually identified clusters from this view, one
can flexibly select the corresponding data entities with lasso-selections
(Figure EmbeddingTree: Hierarchical Exploration of Entity Features in
Embedding-b1) and investigate detailed entity features in the Data Table view.

The Data Table view (Figure EmbeddingTree: Hierarchical Exploration of Entity
Features in Embeddingc) shows the raw data entities in a 2D table, with rows
representing individual entities and columns representing individual entity
features and embedding dimensions. From the table header, users can
conveniently search or filter data entities. They can also sort them by a
given feature of interest.

##  5 Case Studies

This section describes two case studies with embeddings trained from credit-
card transactions and the public 30Music data to show the efficacy of our
EmbeddingTree algorithm and the visualizations.

Figure 2: EmbeddingTree for the embedding of the 30Music data. Entities of the
selected node show a single cluster with random distribution.

###  5.1 Merchant Embedding from Transaction Data

Merchant embedding plays a significant role in many payment-related
applications, such as restaurant recommendation [31] and merchant
identification [32]. The embedding dataset in this study is generated using a
real-world transaction dataset involving $70$ million merchants and $260$
million customers. The merchant embedding is generated by Word2vec [14, 7],
where each merchant is treated as a word and each customer as a document. For
this experiment, we focus on the embeddings for merchants from four cities
(city=0, 1, 2, 3): San Francisco(city=0), Los Angeles (city=1), Chicago
(city=2), and Manhattan (city=3). We choose the three most frequent merchant
category codes (MCC) to filter the embeddings, which are fast food (MCC=5814),
gas station (MCC=5411), and grocery (MCC=5499). After filtering, the numbers
of merchants in the four cities are $9926$, $21456$, $19958$ and $21202$,
respectively. Besides location and MCC, the third feature we are given is the
number of transactions of the merchant (frequency). Since it is not a
categorical feature, we’ve binned it into buckets based on the percentiles
(below $33\%$, $33\%{\sim}66\%$, and above $66\%$). The merchant EmbeddingTree
is built with Algorithm 1 using the splitting criteria introduced in Section
A.1. The resulting tree is shown in Figure EmbeddingTree: Hierarchical
Exploration of Entity Features in Embeddinga.

From the visualization in the Tree view, we can conclude that location is the
most important feature since the tree uses it repetitively in early levels to
split the embedding data. However, under the same city, the features that
further distinguish the entities/merchants are different. For example, when
city=1 (the first left branch in Figure EmbeddingTree: Hierarchical
Exploration of Entity Features in Embeddinga), MCC is the dominant feature,
but it is frequency when city=2 (the corresponding node was collapsed in
Figure EmbeddingTree: Hierarchical Exploration of Entity Features in
Embeddinga). This observation demonstrates the strength of the hierarchical
exploration compared to the parallel exploration. The proposed hierarchical
exploration allows a localized feature importance analysis, and such analysis
is not possible with parallel exploration.

We can also observe the data inconsistency issue from our visualization tool.
The selected leaf node contains $471$ merchants, representing MCC=5499,
City=1, and frequency between $73$ and $101.5$. The embedding should be well
clustered if using only the merchants’ features to do the embedding. However,
we can clearly see the two clusters in Figure EmbeddingTree: Hierarchical
Exploration of Entity Features in Embeddingb. By analyzing the details of the
two clusters, it is interesting to find that most of the merchants on the left
cluster are Square merchants, which accept and process payment with Square.
The merchants on the right are non-Square merchants. The visualization deepens
users’ understanding of the embeddings and helps to discover semantics hidden
in the embeddings.

Additionally, by deploying the visualization tool to be used by domain
practitioners, they notice that our tool can help embedding generation for
unseen entities. Specifically, most embeddings are pre-computed and their
generation relies on the rich interactions between a large number of data
entities. When using the embeddings for downstream tasks, new entities (that
were not seen in the pre-computation stage) may occur, causing the failure of
the pre-trained embeddings. Our EmbeddingTree can help to handle this case.
For example, a newly opened merchant only has a few transactions. Generating
embeddings for such a merchant is very difficult due to the lack of historical
information and interactions with others. By traversing the EmbeddingTree
built for the pre-computed merchant embeddings, we can quickly locate the
embedding cluster that the new merchant should fall into, using its feature
values (e.g. merchant type and location). The mean of all embeddings in the
cluster can be considered as the initial embedding of the new merchant. In
short, EmbeddingTree helps the new merchant find the group it is most similar
to based on its feature information.

###  5.2 User/Track Embedding for Music Recommendation

The embedding examined in this section is learned from the 30Music dataset
[33]. The dataset consists of the listening behavior of $45{,}167$ users and
their meta information, such as gender, age, country, play count, playlist
count, and subscribe type. The dataset also has $5{,}675{,}143$ tracks of
music, and the total number of user-track pairs is $31{,}351{,}945$. We
leverage the listening relationship, i.e., the user-track pairs, to learn the
embedding for individual users using the framework proposed in [26].

The EmbeddingTree built for the resulting users’ embedding is shown in Figure
2. From the hierarchy, the gender of a user is the most significant feature
differentiating their listening behaviors (encoded in the embedding).
Following the top-down exploration, the EmbeddingTree is further split by
playlist, followed by play count.

Investigating users from one leaf node (highlighted in the red circle in
Figure 2), we found that the users’ embedding shows a single cluster in the
Dimensionality Reduction view with no obvious outlier. However, as shown in
the Data Table view, the country and subscribe type of the users have
different values. This reveals a potential deficiency of the learned
embeddings, i.e., they fail to capture the variance from these two features.
To verify this, one can check if the corresponding feature values are correct
or not. For example, the user with id $987$ has a premium subscribe type,
which differs from all other users that have the base type. If the data source
is available, a quick check of this user can help to perform data quality
control. If there is no data quality issue, the subscribe type is probably not
properly encoded in the embedding. The country feature can be explored
similarly. The exploration here helps to disclose potential issues of the
learned embeddings, which motivates users to explicitly inject less-encoded
features into the embedding learning process to refine the embedding.

##  6 Conclusion and Future Works

In this paper, we introduce EmbeddingTree, an approach that allows users to
hierarchically organize embeddings according to their entity features. The
approach helps to analyze the importance of entity features and understand the
relationship between features in the learned embeddings. A visualization tool
is developed and tailored for the visualization and exploration of
hierarchical embeddings. We demonstrate the utility of the visualization tool
with two use cases. One discovers minority features embedded in the merchant
embedding data, the other discloses the feature and embedding inconsistency in
the embedding of the 30Music data. Our algorithm and the visualization tool
have been used by domain practitioners to help their daily practices. In the
future, we would like to conduct statistical user studies to comprehensively
validate the effectiveness of our system. Also, we plan to enable users to
inject their domain knowledge of different data features into the
EmbeddingTree construction process by manually growing/pruning the tree.
Lastly, we will work with more embedding data to further optimize the
generalizability and scalability of our algorithm and visualizations.

0 Research

##  Appendix A Appendix

###  A.1 Details of Section 3.2: 2-GMM Splitting with Embedding-BIC.

One critical component of the Algorithm in Section 3.1 is the criterion used
to select the best splitting feature. The criterion is computed based on the
Fast-BIC for GMMs inspired by [29].

We assume the embedding after PCA projection can be modeled as two mixture
Gaussians. The expectation-maximization (EM) algorithm is used to jointly
estimate all the parameters and latent variables. The latent variables,
$z_{i,j}$, denote the probability that sample $i$ is in cluster $j$. With $N$
as the number of observations and $J$ be the number of Gaussian clusters (in
this case, $J=2$), $z=\\{z_{1,1},z_{1,2},...,z_{N,J-1},z_{N,J}\\}$, the
complete likelihood (including the latent variables) is:

$P(x,\mu,\sigma,w,z)=\prod_{i=1}^{N}\prod_{j=1}^{J}\\{w_{j}\mathcal{N}(x_{n};\mu_{j},\sigma_{j}^{2})\\}^{z_{i,j}}$
(3)

To find the best binary feature that splits the embedding and forms the best
GMM, we go through every feature. Each candidate binary feature splits the
embeddings into two clusters, each is then formulated as a Gaussian. For each
feature, suppose the first $s$ embeddings have feature value $F^{k}=0$ and the
rest $N-s$ embeddings have feature value $F^{k}=1$. We estimate the weights,
means, and variances for both clusters using the maximum likelihood estimation
(MLE).

$\hat{\mu_{1}}=\frac{1}{s}\sum_{i=1}^{s}x_{i},\hskip
14.22636pt\hat{\sigma_{1}}=\frac{1}{s}\sum_{i=1}^{s}||x_{i}-\hat{\mu_{1}}||^{2},\hskip
14.22636pt\hat{w_{1}}=\frac{s}{N},$
$\hat{\mu_{2}}=\frac{1}{N-s}\sum_{i=s+1}^{N}x_{i},\hskip
14.22636pt\hat{\sigma_{2}}=\frac{1}{N-s}\sum_{i=s+1}^{N}||x_{i}-\hat{\mu_{2}}||^{2},\hskip
14.22636pt\hat{w_{2}}=\frac{N-s}{N}.$

In other words, rather than the soft clustering of GMM, our algorithm performs
a hard clustering. Thus, if $x_{i}$ is in cluster $j$, then $z_{i,j}=1$ and
$z_{i,j^{\prime}}=0$ for all $j\neq j^{\prime}$. Given this approximation, the
likelihood can be obtained by summing over the $z$:

$P(x,\mu,\sigma,w)=\sum_{z}\prod_{i=1}^{N}\prod_{j=1}^{J}\\{w_{j}\mathcal{N}(x_{n};\mu_{j},\sigma_{j}^{2})\\}^{z_{i,j}}$
(4)

Note that $z_{(i\in(0,s],j=1)}=z_{(i\in[s+1,N),j=2)}=1$ and $z_{i,j}=0$,
otherwise, the above equation can be simplified to:

$P(x,\mu,\sigma,w)=\prod_{i=1}^{s}w_{1}\mathcal{N}(x_{i};\mu_{1},{\sigma_{1}}^{2})\prod_{i=s+1}^{N}w_{2}\mathcal{N}(x_{i};\mu_{2},{\sigma_{2}}^{2}).$
(5)

Plugging in the MLE for all the parameters, the maximum log-likelihood
function $\hat{L}=\log P(x,\hat{\mu},\hat{\sigma},\hat{w})$ is:

$\hat{L}=\sum_{i=1}^{s}[\log\hat{w_{1}}+\log\mathcal{N}(x_{i};\hat{\mu_{1}},\hat{\sigma_{1}}^{2})]\\\
+\sum_{i=s+1}^{N}[\log\hat{w}_{2}+\log\mathcal{N}(x_{i};\hat{\mu_{2}},\hat{\sigma_{2}}^{2})].$
(6)

Substituting the probability density function of Gaussian
$\mathcal{N}=\frac{1}{\sqrt{2\pi}\sigma}\exp^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}$
into the above equation and simplifying, we get the following expression for
the log-likelihood for any given $s$:

$\hat{L}_{s}=-\frac{s}{2}\log 2\pi\hat{\sigma_{1}}^{2}-\frac{N-s}{2}\log
2\pi\hat{\sigma_{2}}^{2}\\\ +s\log\hat{w_{1}}+(N-s)\log\hat{w_{2}},$ (7)

in which we have dropped terms that are not functions of the parameters. We
call our algorithm Embedding-BIC algorithm. Specifically, it chooses the
feature as the split-point that maximizes $\hat{L_{s}}$. The detailed
pseudocode for this algorithm can be found in Algorithm 1.

Algorithm 2 Calculate the BIC score for embedding $\mathbf{X}$ with a given
binary feature $\mathbf{F}^{k}$.

1:procedure Embedding-BIC($\mathbf{X},\mathbf{F}^{k})$)

2: for $i\in[1,N]$ do

3: $x_{i}=PCA(\mathbf{x}_{i})$

4: $C_{1}\leftarrow\mathbf{f}_{i}^{k}=0$

5: $C_{2}\leftarrow\mathbf{f}_{i}^{k}=1$

6: for $j\in\\{1,2\\}$ do

7: $n_{j}=length(C_{j})$

8: $\hat{w_{j}}=n_{j}/n$

9: $\hat{\mu}_{j}=\frac{1}{n_{j}}\sum_{x_{i}\in C_{j}}x_{i}$

10: $\hat{\sigma_{j}}^{2}=\frac{1}{n_{j}}\sum_{x_{i}\in
C_{j}}(x_{i}-\hat{\mu_{j}})^{2}$

11: BIC_curr $=-\frac{n_{1}}{2}\log
2\pi\hat{\sigma}_{1}^{2}-\frac{n_{2}}{2}\log
2\pi\hat{\sigma}_{2}^{2}+n_{1}\log\hat{w}_{1}+n_{2}\log\hat{w}_{1}$

12: return BIC_curr

## References

  * [1] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” _arXiv:1301.3781_ , 2013.
  * [2] J. Camacho-Collados and M. T. Pilehvar, “From word to sense embeddings: A survey on vector representations of meaning,” _JAIR_ , vol. 63, pp. 743–788, 2018.
  * [3] V. Tshitoyan, J. Dagdelen, L. Weston, A. Dunn, Z. Rong, O. Kononova, K. A. Persson, G. Ceder, and A. Jain, “Unsupervised word embeddings capture latent knowledge from materials science literature,” _Nature_ , vol. 571, no. 7763, p. 95, 2019.
  * [4] M. Grbovic and H. Cheng, “Real-time personalization using embeddings for search ranking at airbnb,” in _KDD_ , 2018, pp. 311–320.
  * [5] J. Wang, P. Huang, H. Zhao, Z. Zhang, B. Zhao, and D. L. Lee, “Billion-scale commodity embedding for e-commerce recommendation in alibaba,” in _KDD_ , 2018, pp. 839–848.
  * [6] X. Zhao, R. Louca, D. Hu, and L. Hong, “Learning item-interaction embeddings for user recommendations,” _arXiv:1812.04407_ , 2018.
  * [7] M. Du, R. Christensen, W. Zhang, and F. Li, “Pcard: Personalized restaurants recommendation from card payment transaction records,” in _WWW_ , 2019, pp. 2687–2693.
  * [8] D. Smilkov, N. Thorat, C. Nicholson, E. Reif, F. B. Viégas, and M. Wattenberg, “Embedding projector: Interactive visualization and interpretation of embeddings,” _arXiv preprint arXiv:1611.05469_ , 2016.
  * [9] F. Heimerl, C. Kralj, T. Moller, and M. Gleicher, “embcomp: Visual interactive comparison of vector embeddings,” _IEEE Transactions on Visualization and Computer Graphics_ , 2020.
  * [10] Y. Liu, E. Jun, Q. Li, and J. Heer, “Latent space cartography: Visual analysis of vector space embeddings,” in _Computer Graphics Forum_ , vol. 38, no. 3. Wiley Online Library, 2019, pp. 67–78.
  * [11] S. Liu, P.-T. Bremer, J. J. Thiagarajan, V. Srikumar, B. Wang, Y. Livnat, and V. Pascucci, “Visual exploration of semantic relationships in neural word embeddings,” _IEEE transactions on visualization and computer graphics_ , vol. 24, no. 1, pp. 553–562, 2017.
  * [12] M. Berger, “Visually analyzing contextualized embeddings,” in _2020 IEEE Visualization Conference (VIS)_. IEEE, 2020, pp. 276–280.
  * [13] Y. Wang, Y. Zheng, Y. Peng, M. Yeh, Z. Zhuang, D. Mahashweta, B. Mangesh, F. Li, W. Zhang, and J. M. Phillips, “Constrained non-affine alignment of embeddings,” in _2021 IEEE International Conference on Data Mining (ICDM)_. IEEE, 2021, pp. 1403–1408.
  * [14] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” vol. 26, 2013\.
  * [15] J. Pennington, R. Socher, and C. D. Manning, “Glove: Global vectors for word representation,” in _Empirical Methods in Natural Language Processing (EMNLP)_ , 2014, pp. 1532–1543.
  * [16] Z. Harris, “Distributional structure,” _Word_ , vol. 10, no. 23, pp. 146–162, 1954.
  * [17] B. Li and D. Pi, “Learning deep neural networks for node classification,” _Expert Systems with Applications_ , vol. 137, pp. 324–334, 2019.
  * [18] W. Fan, Y. Ma, Q. Li, Y. He, E. Zhao, J. Tang, and D. Yin, “Graph neural networks for social recommendation,” in _WWW_ , 2019, pp. 417–426.
  * [19] X. Wei, L. Xu, B. Cao, and P. S. Yu, “Cross view link prediction by learning noise-resilient representation consensus,” in _WWW_ , 2017, pp. 1611–1619.
  * [20] H. Cai, V. W. Zheng, and K. C.-C. Chang, “A comprehensive survey of graph embedding: Problems, techniques, and applications,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 30, no. 9, pp. 1616–1637, 2018.
  * [21] P. Goyal and E. Ferrara, “Graph embedding techniques, applications, and performance: A survey,” _Knowledge-Based Systems_ , vol. 151, pp. 78–94, 2018.
  * [22] T. Rahman, B. Surma, M. Backes, and Y. Zhang, “Fairwalk: Towards fair graph embedding,” in _Proceedings of the 2019 International Joint Conferences on Artifical Intelligence (IJCAI)_ , 2019.
  * [23] Q. Wang, Z. Mao, B. Wang, and L. Guo, “Knowledge graph embedding: A survey of approaches and applications,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 29, no. 12, pp. 2724–2743, 2017.
  * [24] W. Hamilton, Z. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [25] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” in _Advances in neural information processing systems_ , 2013, pp. 3111–3119.
  * [26] C.-C. M. Yeh, D. Gelda, Z. Zhuang, Y. Zheng, L. Gou, and W. Zhang, “Towards a flexible embedding learning framework,” in _2020 International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2020, pp. 605–612.
  * [27] X. Liu and J. Wang, “Latentvis: Investigating and comparing variational auto-encoders via their latent space.” in _CIKM Workshops_ , 2020.
  * [28] Q. Li, K. S. Njotoprawiro, H. Haleem, Q. Chen, C. Yi, and X. Ma, “Embeddingvis: A visual analytics approach to comparative network embedding inspection,” in _2018 IEEE Conference on Visual Analytics Science and Technology (VAST)_. IEEE, 2018, pp. 48–59.
  * [29] A. Boggust, B. Carter, and A. Satyanarayan, “Embedding comparator: Visualizing differences in global structure and local neighborhoods via small multiples,” in _27th international conference on intelligent user interfaces_ , 2022, pp. 746–766.
  * [30] M. Madhyastha, G. Li, V. Strnadová-Neeley, J. Browne, J. T. Vogelstein, R. Burns, and C. E. Priebe, “Geodesic forests,” in _KDD_ , 2020, pp. 513–523.
  * [31] M. Du, R. Christensen, W. Zhang, and F. Li, “Pcard: Personalized restaurants recommendation from card payment transaction records,” in _WWW_ , 2019, pp. 2687–2693.
  * [32] C.-C. M. Yeh, Z. Zhuang, Y. Zheng, L. Wang, J. Wang, and W. Zhang, “Merchant category identification using credit card transactions,” in _IEEE International Conference on Big Data_ , 2020, pp. 1736–1744.
  * [33] R. Turrin, M. Quadrana, A. Condorelli, R. Pagano, and P. Cremonesi, “30music listening and playlists dataset.” in _RecSys Posters_ , 2015.