# Charge Singlets and Orbital Selective Charge Density Wave Transitions

Yuxi Zhang  Department of Physics, University of California, Davis, CA
95616,USA  Chunhan Feng  Department of Physics, University of California,
Davis, CA 95616,USA  Rubem Mondaini  Beijing Computational Science Research
Center, Beijing 100084, China  G.G. Batrouni  Université Côte d’Azur, CNRS,
Institut de Physique de Nice (INPHYNI), 06000 Nice, France  Centre for Quantum
Technologies, National University of Singapore, 2 Science Drive 3, 117542
Singapore  Department of Physics, National University of Singapore, 2 Science
Drive 3, 117542 Singapore  Beijing Computational Science Research Center,
Beijing 100084, China  Richard T. Scalettar  Department of Physics, University
of California, Davis, CA 95616,USA

###### Abstract

The possibility of ‘orbitally selective Mott transitions’ within a multi-band
Hubbard model, in which one orbital with large on-site electron-electron
repulsion $U_{1}$ is insulating and another orbital, to which it is
hybridized, with small $U_{-1}$, is metallic, is a problem of long-standing
debate and investigation. In this paper we study an analogous phenomenon, the
co-existence of metallic and insulating bands in a system of orbitals with
different electron-phonon coupling. To this end, we examine two variants of
the bilayer Holstein model: a uniform bilayer and a ‘Holstein-metal interface’
where the electron-phonon coupling (EPC), $\lambda$, is zero in the ‘metallic’
layer. In the uniform bilayer Holstein model, charge density wave (CDW) order
dominates at small interlayer hybridization $t_{3}$, but decreases and
eventually vanishes as $t_{3}$ grows, providing a charge analog of singlet
(spin liquid) physics. In the interface case, we show that CDW order
penetrates into the metal layer and forms long-range CDW order at intermediate
ratio of inter- to intra-layer hopping strengths, $1.4\lesssim t_{3}/t\lesssim
3.4$. This is consistent with the occurrence of an ‘orbitally selective CDW’
regime at weak $t_{3}$ in which the layer with $\lambda_{1}\neq 0$ exhibits
long range charge order, but the ‘metallic layer’ with $\lambda_{-1}=0$, to
which it is hybridized, does not.

##  I Introduction

Over the last several decades, much attention has been focused in the
condensed matter community on layered materials. One prominent example is that
of the cuprate superconductors (SC) [1, 2, 3, 4, 5]. Bilayer graphene [6, 7,
8, 9, 10, 11] is another, more recent, realization. From a theoretical
perspective, bilayer materials offer an opportunity to explore the competition
between the formation of long range order at weak interlayer coupling and
collections of independent local degrees of freedom in the limit of strong
interlayer coupling. Computational studies have lent considerable insight into
these phenomena, including quantitative values for the quantum critical points
[12, 13, 14, 15, 16] separating antiferromagnetic and singlet phases at zero
temperature.

This competition is central to that which occurs in multiorbital systems,
notably the interplay of Ruderman-Kittel-Kasuya-Yosida order and singlet
formation in the Kondo lattice and periodic Anderson models [17, 18, 19]. This
close analogy originates in the observation that, in calculations on a model
Hamiltonian, there is no difference between multi-layer and multi-orbital
descriptions, apart from the interpretation of the additional label of the
fermionic species. For this reason we will use the two terminologies
interchangeably here. In multi-orbital language, one of the key conceptual
interests is the possibility that the distinct values of the ratio of
interaction strength to kinetic energy in the different bands might result in
separate insulator transitions, i.e. the possibility of an ‘orbital-selective’
Mott transition (OSMT) [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

Here we study analogous questions concerning bilayer (bi-orbital) systems in
which the fermions interact with phonon degrees of freedom rather than via
direct electron-electron correlations. A precise mathematical description of
the mapping between the two situations is discussed. Quantum Monte Carlo (QMC)
simulations have already been applied to the analysis of charge-density wave
(CDW) and superconducting (SC) transitions in the single band Holstein model
[32]. However, thus far, work has focused mostly on two-dimensional or three-
dimensional models with a single kinetic energy scale [33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46].

Using QMC simulations of the two-band Holstein model at half-filling, we will
address the following questions concerning the effects of interband
hybridization $t_{3}$: (i) Is there a transition in which CDW order is
destroyed as $t_{3}$ is increased? What is the value of the critical coupling
associated with the quantum critical point (QCP) in the ground state and the
critical temperature for the thermal transitions at finite $T$? (ii) In a
situation where the electron-phonon energy scales in the two bands are very
different, can CDW order in one band coexist with metallic behavior in the
other? These issues are in direct analogy with those addressed in multiband
Hubbard Hamiltonians; we will discuss similarities and differences between the
resulting phenomena.

The paper is organized as follows: Sec. II starts by describing the model
Hamiltonian, with Sec. III introducing the numerical algorithm employed. The
main results of both the bilayer Holstein model and the Holstein-metal
interface are given in Sec. IV and Sec. V, respectively. A contrast between
the two models is drawn in Sec. VI, by analyzing local quantities and gaps to
excitations; Sec. VII summarizes our results.

##  II Layered Holstein Hamiltonian

We focus on the bilayer Holstein model

$\displaystyle\mathcal{\hat{H}}=$ $\displaystyle-\sum_{\langle
ij\rangle,l,\sigma}\big{(}t^{\phantom{\dagger}}_{l}\,\hat{c}^{\dagger}_{il\sigma}\hat{c}^{\phantom{\dagger}}_{jl\sigma}+{\rm
h.c.}\big{)}-\sum_{i,l,\sigma}\mu_{l}\hat{n}_{il\sigma}$
$\displaystyle+\frac{1}{2M}\sum_{il}\hat{P}^{2}_{il}+\frac{1}{2}\sum_{i,l}\omega_{l}^{2}\hat{X}^{2}_{il}+\sum_{i,l,\sigma}\lambda_{l}\,\hat{n}_{il\sigma}\hat{X}_{il}$
$\displaystyle-\sum_{i,\langle
ll^{\prime}\rangle,\sigma}\big{(}t^{\phantom{\dagger}}_{ll^{\prime}}\,\hat{c}^{\dagger}_{il\sigma}\hat{c}^{\phantom{\dagger}}_{il^{\prime}\sigma}+{\rm
h.c.}\big{)}\,\,.$ (1)

$\hat{c}^{\phantom{\dagger}}_{il\sigma}(\hat{c}^{\dagger}_{il\sigma})$ are
annihilation (creation) operators for an electron on layer $l(=\pm 1)$, site
$i$ with spin $\sigma$, and
$\hat{n}_{il\sigma}=\hat{c}^{\dagger}_{il\sigma}\hat{c}^{\phantom{\dagger}}_{il\sigma}$
is the number operator. $t_{l}$ and $t_{ll^{\prime}}=t_{3}$ denote the intra-
and inter-layer hopping respectively. Phonons are represented by local
(dispersionless) quantum harmonic oscillators with frequency $\omega_{l}$, and
on-site electron-phonon interaction on layer $l$ is introduced via
$\lambda_{l}$. We choose intralayer hopping $t_{l}=t=1$ throughout this work
to set the energy scale, and all simulations are done at half-filling
$\langle\hat{n}_{il}\rangle=1$, which can be achieved by setting the chemical
potential $\mu_{l}=-\lambda_{l}^{2}/\omega_{l}^{2}$; phonon mass is set as
$M=1$. Each layer is an $L\times L$ site square lattice, as sketched in the
inset of Fig. 1(b), with $N=2\times L\times L$ being the total number of
sites. We focus on two cases in this work: a uniform bilayer Holstein model
where $t_{1}=t_{-1}=t$, $\mu_{1}=\mu_{-1}=\mu$,
$\omega_{1}=\omega_{-1}=\omega$ and $\lambda_{+1}=\lambda_{-1}=\lambda$; and
an interface between Holstein layer and “metal” layer, where only layer $l=+1$
has a non-zero electron-phonon coupling $\lambda_{+1}\neq 0$ and layer $l=-1$
has $\lambda_{-1}=0$. We employ a recently developed Langevin quantum Monte
Carlo (QMC) method [47] discussed in the next section.

We first define the local observables including the (layer-dependent) double
occupancy,

$\displaystyle{\cal
D}_{l}\equiv\langle\hat{n}_{il\uparrow}\hat{n}_{il\downarrow}\rangle$ (2)

the near-neighbor intra-layer Green’s function,

$\displaystyle{\cal G}_{\langle ij\rangle
l}\equiv-\langle\hat{c}^{\dagger}_{il\sigma}\hat{c}^{\phantom{\dagger}}_{jl\sigma}+\hat{c}^{\dagger}_{jl\sigma}\hat{c}^{\phantom{\dagger}}_{il\sigma}\rangle\,,$
(3)

and the near-neighbor inter-layer Green’s function,

$\displaystyle{\cal G}_{\langle
ll^{\prime}\rangle}\equiv-\langle\hat{c}^{\dagger}_{il\sigma}\hat{c}^{\phantom{\dagger}}_{il^{\prime}\sigma}+\hat{c}^{\dagger}_{il^{\prime}\sigma}\hat{c}^{\phantom{\dagger}}_{il\sigma}\rangle\,.$
(4)

When multiplied by their associated hopping integrals, $t\,{\cal G}_{\langle
ij\rangle l}$ and $t_{3}\,{\cal G}_{\langle ll^{\prime}\rangle}$ give the
intra- and inter-layer kinetic energies per site.

Two further observables, the density-density and pair-pair correlators, aid in
characterizing the excitations between the planes.

$\displaystyle d_{-1,1}$
$\displaystyle\equiv\frac{1}{4}\langle\hat{n}^{\phantom{\dagger}}_{i,1}\hat{n}^{\phantom{\dagger}}_{i,-1}-1\rangle$
$\displaystyle p_{-1,1}$
$\displaystyle\equiv-\frac{1}{4}\langle\hat{\Delta}^{\phantom{\dagger}}_{i,1}\hat{\Delta}^{\dagger}_{i,-1}+\hat{\Delta}^{\dagger}_{i,1}\hat{\Delta}^{\phantom{\dagger}}_{i,-1}\rangle$
$\displaystyle\hat{\Delta}^{\dagger}_{il}$
$\displaystyle\equiv\hat{c}^{\dagger}_{il\uparrow}\hat{c}^{\dagger}_{il\downarrow}.$
(5)

$d_{-1,1}$ and $p_{-1,1}$ are the analogs of the $zz$ and $xy$ spin
correlations which enter into the characterization of interlayer singlet
formation in the Hubbard and Heisenberg bilayers, see Appendix A. Because of
rotational symmetry of those models, their magnetic analogs, obtained by the
transformation
$\hat{c}^{\phantom{\dagger}}_{il\downarrow}\rightarrow\hat{c}^{\dagger}_{il\downarrow}$
are identical in value. $d_{-1,1}=p_{-1,1}$ would also hold in the attractive
Hubbard Hamiltonian. Here, in the Holstein model, rotational symmetry is
broken and we have $d_{-1,1}\neq p_{-1,1}$. We will discuss the implications
further in the sections to follow.

Characterization of the CDW formation in the thermodynamic limit can be made
by the analysis of the (layer-resolved) structure factor,

$\displaystyle S^{\rm
cdw}_{l}=\frac{2}{N}\sum_{ij}(-1)^{i+j}\langle\hat{n}_{il}\hat{n}_{jl}\rangle,$
(6)

with $\hat{n}_{il}=\sum_{\sigma}\hat{n}_{il\sigma}$. $S^{\rm cdw}_{l}$ samples
correlations across the entire lattice, and hence is a primary tool in the
determination of long range order.

In the case of the uniform bilayer, the quantities defined in Eqs. 2-4 and 6
are independent of the layer index $l$, and in this case we suppress this
index. But for the ‘interface’ geometry, which includes one layer with
$\lambda_{1}\neq 0$ and another with $\lambda_{-1}=0$, measurements performed
on the two layers are inequivalent.

The layer-resolved single-particle spectral function $A_{l}(\omega)$ is
obtained by using the maximum entropy method to invert the integral equation
relating the imaginary time dependent Green’s function $G_{i=0}(\tau)$ and
$A(\omega)$:

$\displaystyle G_{i=0}(\tau)$ $\displaystyle=\int
d\omega\frac{e^{-\tau\omega}}{1+e^{\beta\omega}}\,A(\omega)$ $\displaystyle
G_{i}(\tau)$
$\displaystyle=\langle\hat{c}^{\phantom{\dagger}}_{i}(\tau)\hat{c}^{\dagger}_{0}(0)\rangle=\langle
e^{\tau\mathcal{H}}\hat{c}^{\phantom{\dagger}}_{i}(0)e^{-\tau\mathcal{H}}\hat{c}^{\dagger}_{0}(0)\rangle\,.$
(7)

$\tau$ represents imaginary time; layer and spin indices are omitted here for
simplicity. The appropriate local $G$ is used to get $A_{l}(\omega)$ for; each
layer $l$ in the interface geometry.

We advance our key results in Figures 1(a) and 1(b): (i) At weak $t_{3}$ there
is a phase transition at finite temperature $T_{c}$ to a state with long range
charge order. In the bilayer case, $T_{c}$ initially increases with $t_{3}$ as
the charge order is enhanced by increased coordination number. (ii) At $T=0$,
in both the Holstein bilayer and the Holstein-metal interface, CDW order is
destroyed for $t_{3}$ exceeding a quantum critical value. (iii) The phase
diagram of the interface geometry exhibits an ‘orbitally selective CDW phase’
(OSCDW) at low $T$ and weak $t_{3}$. The specific description of how these
phase diagrams are obtained is given in the corresponding section containing
the main results of each model.

Figure 1:  (a) Phase diagram of the Holstein bilayer giving the CDW transition
temperature $T_{c}$ as a function of inter-layer hopping $t_{3}$. Two values
of electron-phonon coupling, $\lambda=0.2$ and $\lambda=0.16$, are shown. Data
points are obtained by finite size scaling (FSS) analysis as shown in Fig. 3.
Dashed lines are guides to the eye. Inset shows crossing plot of $S^{\rm
cdw}/L^{7/4}$ versus $t_{3}$ at $\lambda=0.2$ and low temperature $\beta=20$
and determines the quantum critical value for the CDW-charge singlet
transition. (b) Analog of (a), but for the Holstein-metal interface. The four
data points (blue crosses) at $T=0.1$ and $T=0.05$ are obtained by analysis of
the scaled structure factor, Fig. 6. Inset shows a sketch of a bilayer with
relevant terms in Eq. (II) marked. QCPs are marked by stars on the
$t_{3}$-axes in panels (a) and (b). The phonon frequency $\omega=0.1$ for all
data.

##  III Langevin Simulation Algorithm

We employ a recently developed Fourier accelerated Langevin quantum Monte
Carlo (QMC) method [47]. The partition function of the Holstein Hamiltonian is
written as a path integral $\mathcal{Z}=\Tr e^{-\beta\mathcal{\hat{H}}}=\Tr
e^{-\Delta\tau\mathcal{\hat{H}}}e^{-\Delta\tau\mathcal{\hat{H}}}\cdots
e^{-\Delta\tau\mathcal{\hat{H}}}$ where the inverse temperature
$\beta=L_{\tau}\Delta\tau$ is discretized along the “imaginary time” axis.
Complete sets of phonon eigenstates $|\\{x_{i,\tau}\\}\rangle$ are inserted at
each time slice, allowing the action of the phonon operators to be evaluated.
In so doing, we convert the quantum problem into a classical problem in one
higher dimension. Since the fermion operators appear only as quadratic forms,
we can trace over the associated degrees of freedom, leaving the partition
function dependent only on the phonon field $\\{x_{i,\tau}\\}$,

$\displaystyle\mathcal{Z}=\int\mathcal{D}x_{i,\tau}e^{-S_{\rm ph}}[\det
M(\\{x_{i,\tau}\\})]^{2}=\int\mathcal{D}x_{i,\tau}e^{-S}\,,$ (8)

where the “phonon action”

$\displaystyle S_{\rm
ph}=\frac{\Delta\tau}{2}\left[\omega^{2}\sum_{i}x_{i,\tau}^{2}+\sum_{i}\left(\frac{x_{i,\tau+1}-x_{i,\tau}}{\Delta\tau}\right)^{2}\right]$
(9)

and

$\displaystyle S=S_{\rm ph}-\ln(\det M)^{2}\,.$ (10)

Here $M$ is a sparse matrix of dimension $NL_{\tau}$ whose detailed form is
given in Ref. [47]. The square of the determinant appears because up and down
fermionic species have the same coupling to the phonons. As a consequence,
there is no sign problem [48, 49]. In order to sample the phonon coordinates,
instead of using the usual Metropolis algorithm, we evolve $\\{x_{i,\tau}\\}$
using the discretized Langevin equation, whose simplest form is given by the
first order Euler discretization,

$\displaystyle x_{i,\tau,t+dt}=x_{i,\tau,t}-dt\frac{\partial S}{\partial
x_{i,\tau,t}}+\sqrt{2\,dt}\,\eta_{i,\tau,t}\,,$ (11)

where t is the Langevin time, and $\eta$ is a Gaussian distributed stochastic
variable. In practice, in our simulations we make use of a higher order Runge-
Kutta discretization [47] which reduces the discretization error to ${\cal
O}(dt^{2})$. Throughout this work, the Langevin time step dt is chosen as
0.002, which has been shown to be sufficiently small so that the Langevin time
discretization error is of the same order, or smaller than, statistical errors
in typical simulations [47]. It can be demonstrated that, in the stationary
limit, this Markov process generates configurations which are drawn from the
exponential of the action of Eq. 10. The computational kernel is the
calculation of the partial derivatives of the action via

$\displaystyle\frac{\partial S}{\partial x_{i,\tau,t}}=\frac{\partial S_{\rm
ph}}{\partial x_{i,\tau,t}}-2\Tr{\frac{\partial M}{\partial
x_{i,\tau,t}}M^{-1}},$ (12)

where the trace is evaluated using a stochastic estimator [47]. Comparing to
the conventional determinant quantum Monte Carlo (DQMC) method, which is an
$O(N^{3}L_{\tau})$ approach, the Langevin method scales as $O(NL_{\tau})$
(although with a larger prefactor, so that there is a cross-over $N$ at which
the Langevin approach becomes the more efficient method). This enables
simulations to reach considerably larger lattice sizes. In this work, we
analyze systems up to $N=800$ sites. The efficiency of the Langevin approach
results from the sparsity of the matrix $M$ and the fact that computing the
action of $M^{-1}$ on a vector can be done iteratively in a number of steps
which does not grow with $N$ [47], with appropriate pre-conditioning. The
Langevin dynamics we employ is particularly effective in the adiabatic limit
of small phonon-frequencies, where the density of zeros of individual fermion
determinants is negligible [50]. In what follows we fix $\omega=0.1$, thus
simulations are stable, and statistical convergence is quickly obtained over
the course of Markov generation.

##  IV Holstein bilayer

We initially consider two identical layers with $\lambda=0.2,\,\omega=0.1$ and
the question of the destruction of CDW order via the formation of charge
singlets at large interlayer hopping $t_{3}$ before tackling the more complex
issue of selective CDW transitions.

Figure 2:  (a) Charge structure factor $S^{\rm cdw}$; and (b) $p_{-1,1}$ and
$d_{-1,1}$ as a function of interlayer hopping $t_{3}$ for the Holstein
bilayer at low temperature $\beta=16$ and $\lambda_{+1}=\lambda_{-1}=0.2$.
$S^{\rm cdw}$ shows significant finite size effects in the ordered phase
$t_{3}\lesssim 4.8$. Note that $d_{-1,1}$ vanishes at $t_{3}=0$, but jumps
discontinuously to a non-zero value for infinitesmal $t_{3}$. Panels (c,d) are
analog of (a,b) for the Holstein-metal interface. The two curves in (c)
correspond to layers $l=+1$ and $l=-1$, with $\lambda_{+1}=0.2$ and
$\lambda_{-1}=0$ respectively and temperature $\beta=20$. In all plots the
phonon frequency is set at $\omega=0.1$.

Figure 2(a) gives the CDW structure factor $S^{\rm cdw}$ as a function of
$t_{3}$ at low temperature for two lattice sizes. Below $t_{3,c}\approx 4.8$,
$S^{\rm cdw}$ is large, and grows with lattice size, suggesting long range
charge order. Figure 2(b) focuses on the interlayer density-density $d_{-1,1}$
and pair-pair $p_{-1,1}$ correlations. For small $t_{3}$, only $d_{-1,1}$ is
large in magnitude, indicating coherence in the charge order between the two
layers. As $t_{3}$ increases, intersheet pair correlations $p_{-1,1}$ develop.
$d_{-1,1}$ and $p_{-1,1}$ then become nearly degenerate at $t_{3,c}$,
signalling the loss of CDW order and entry into the ‘charge singlet’ phase.
Together, Figs. 2(a) and 2(b) motivate the bilayer phase diagram of Fig. 1(a).

Although the QCP in this Holstein bilayer is closely analogous to that
occurring in Hubbard and Heisenberg bilayers as well as the periodic Anderson
model, in those cases the electron-electron interaction gives rise to magnetic
phases which form due to the breaking of a continuous spin symmetry. Thus in
2D and quasi-2D geometries, no long range order is possible at finite $T$. In
contrast, here for the Holstein model, charge and pairing order are not
degenerate, as emphasized by the data of Fig. 2(b). CDW correlations dominate
at half-filling and a finite temperature phase transition can occur,
terminating at a QCP as shown in Figs. 1(a) and 1(b). This distinction means
that, in principle, our characterization of the unordered phase as a ‘charge
singlet’ is somewhat loose: in the usual spin singlet the $x,y,z$ components
of the spin-spin correlations on the two layers (or in the two orbitals) are
equal. With that said, the equivalence of $p^{\phantom{\dagger}}_{-1,1}$ and
$d^{\phantom{\dagger}}_{-1,1}$ in the large $t_{3}$ regime points to an
emergent restoration of the symmetry (see Appendix A). It is worth noting that
in the absence of $t_{3}$, e.g. in the 2D Holstein model, this restoration
does not occur until the anti-adiabatic limit is reached, which requires very
large values of $\omega$ [51].

Figure 3: (a) CDW structure factor $S^{\rm cdw}$ dependence on the inverse
temperature $\beta$ and finite size scaling of the Holstein bilayer at
$t_{3}=2$. Both the crossing plot (b), and the full data collapse (c) using 2D
Ising critical exponents and $\beta_{c}\simeq 3.0$.

Figure 3 provides details of the behavior of the CDW structure factor. The top
panel (a) gives raw values for $S^{\rm cdw}$ as a function of $\beta$ at
$t_{3}=2$ for different lattice sizes. At low $\beta$ (high temperature), the
correlation length $\xi$ is short and $S^{\rm cdw}$ is independent of $L$. As
$\beta$ increases, so does $\xi$ and when $\xi\sim L$, $S^{\rm cdw}$ becomes
sensitive to $L$. This separation of the curves provides a crude estimate for
$\beta_{c}$, which may then be determined precisely by finite size scaling
(FSS).

In particular, in the vicinity of the critical temperature $T_{c}$, the CDW
structure factor measured on finite lattices of linear dimension $L$ should
obey,

$\displaystyle S^{\rm cdw}\sim
L^{\gamma/\nu}f\left(\frac{T-T_{c}}{T_{c}}L^{1/\nu}\right)\,\,.$ (13)

As a consequence, when plotting $S^{\rm cdw}/L^{\gamma/\nu}$ as a function of
the inverse temperature $\beta$, different sizes $L$ cross at
$\beta=\beta_{c}$ [Fig. 3 (b)]. Following the scaling form given in Eq. (13)
we note that when plotted against $(\beta-\beta_{c})L^{1/\nu}$ all data
collapse on a single curve – see Figs. 3(c) (and later for the Holstein-metal
interface in Fig. 5). In this analysis we have used the critical exponents of
the 2D Ising universality class ($\gamma=7/4$ and $\nu=1$), since the CDW
phase breaks a $Z_{2}$ symmetry. A discussion of the degree to which the
collapse worsens, and hence the accuracy with which the exponents can be
determined, is given in the Appendix B.

Using such scaling procedure for various values of the interplane
hybridization, allows us to extract the location of the thermal transition, as
compiled in Fig. 1(a), using two values of the electron-phonon interaction. In
this geometry, the critical temperature $T_{c}$ initially increases as a
consequence of the larger coordination number when the planes are coupled -
the 2D to 3D crossover. However, at large $t_{3}$ the critical temperature
decreases and ultimately vanishes at a quantum critical point.

##  V Holstein-Metal interface

We next consider the ‘Holstein-metal interface’ in which layer $l=+1$ has
nonzero $\lambda_{+1}$ but $\lambda_{-1}=0$. The two layers are in contact via
hybridization. Here, in addition to the question of charge singlet formation
at large $t_{3}$, quenching CDW order, a different fundamental question
arises: to what extent do CDW correlations in layer $l=+1$ ‘penetrate’ into
layer $l=-1$, and, conversely, is the CDW in layer $l=+1$ disrupted by contact
with the ‘metallic’ layer? We choose $\lambda=0.2$ and $\omega=0.1$ as in the
previous section.

Figure 2(c) shows the CDW structure factor $S^{\rm cdw}$ in the two layers.
$S^{\rm cdw}_{+1}$ decreases steadily with $t_{3}$: additional quantum
fluctuations associated with contact with the metal reduce charge order. In
contrast, $S^{\rm cdw}_{-1}$ is non-monotonic: charge order is initially
induced in the metal via contact with the Holstein layer, but ultimately large
$t_{3}$ is inimical to it. The behavior of $S^{\rm cdw}_{-1}$ provides a first
clue that order in layer $l=-1$ might occur only for intermediate $t_{3}$.
Figure 2(d) gives the interlayer density-density $d_{-1,1}$ and pair-pair
$p_{-1,1}$ correlations for this interface geometry. The primary difference
from the original bilayer case is the gradual development of $d_{-1,1}$ with
$t_{3}$. This is a consequence of the absence of CDW in the metal layer when
$t_{3}$ vanishes. The interlayer hopping thus must not only couple the charge
correlations, but also induce them in layer $l=-1$. Similar to the bilayer
case, $d_{-1,1}$ and $p_{-1,1}$ become degenerate for large $t_{3}$. This is
again a signature of entering into the charge singlet phase.

Figure 4: (a) CDW structure factor, $S^{\rm cdw}$, dependence on the inverse
temperature, $\beta$, for layer $l=1$ of the Holstein interface at $t_{3}=1$;
(b) using 2D Ising critical exponents for finite size scaling (FSS). Panels
(c) and (d) display the same but for the metallic layer $l=-1$. Panels (e–h)
display the corresponding data for $t_{3}=2$. FSS in (f) and (h) show the same
critical temperature for both layers at $t_{3}=2$, in contrast to $t_{3}=1$,
where layer $l=+1$ (b) exhibits a clear CDW transition whereas data for layer
$l=-1$ (d) does not exhibit crossing when using Ising critical exponents.
Figure 5: Full data collapse of the scaled CDW order parameter versus scaled
reduced temperature in the ‘Holstein interface’ system. Only the Holstein
layer (layer 1) shows a single universal curve for $t_{3}=1$ (a), while both
collapse for $t_{3}=2t$ (panels c,d).

We now turn to a more careful FSS study of the layer-resolved $S^{\rm
cdw}_{l}$. Our main interest is in determining how long range order in the two
layers evolves with $t_{3}$. Figure 4 displays a detailed analysis of two
representative values, $t_{3}=1$ and $t_{3}=2$. The former is a case when
$S^{\rm cdw}_{-1}$ is just beginning to develop, and the latter is when
$S^{\rm cdw}_{-1}$ has reached its maximal induced value [see Fig. 2(c)].
There is a superficial resemblance in the unscaled data for both values of
$t_{3}$, which rise as the temperature is lowered ($\beta$ increases) and also
increase with system size. A proper scaling analysis, however, reveals a
profound distinction. As seen in Figs. 4(b) and 4(f), for both values of
$t_{3}$, the layer $l=+1$ with non-zero electron-phonon coupling
$\lambda_{+1}=0.2$, has a scaled structure factor $L^{-7/4}S^{\rm cdw}_{+1}$
which exhibits a sharp crossing, indicating a finite temperature transition to
long range CDW order. When $t_{3}=2$, this crossing occurs for the metallic
layer with $\lambda_{-1}=0$ as well [Fig. 4(h)]. However, when $t_{3}=1$ the
data for the metallic layer do not cross for the studied system sizes [Fig.
4(d)], namely $L=8\,$-$20$ for both $t_{3}$ values: The $L=12$,$16$ and $20$
curves converge at $\beta>5.4$ instead of crossing. The $L=8$ data do not
scale with the other lattice sizes at all. This distinction becomes even more
apparent in Fig. 5, where a simultaneous data collapse for the scaled
structure factor can be made possible at the same temperature for $t_{3}=2$,
while it is unattainable for $t_{3}=1$.

Figure 6: Crossing plot of $S^{\rm cdw}_{-1}/L^{7/4}$ (metallic layer) versus
$t_{3}$ at $\lambda=0.2$ and $\omega=0.1$, at low temperature (a) $\beta=10$
(b) $\beta=20$, for the Holstein-metal interface. Two crossings are seen for
in each case: (a) $t_{3}\sim 1.3$ and $t_{3}\sim 3.2$ ($\beta=20$) and (b)
$t_{3}\sim 1.33$ and $t_{3}\sim 3.1$ ($\beta=10$). All crossings are obtained
with the critical exponents of the two-dimensional Ising model consistent with
breaking the discrete $Z_{2}$ symmetry. The metallic layer, $l=-1$, is in the
CDW phase only between the two values of $t_{3}$ where the curves cross.

We conclude that for $t_{3}=2$, the interface geometry has CDW order in both
layers, with long range correlations in the metallic layer induced by
proximity to the Holstein layer. For $t_{3}=t$, the interface geometry
exhibits orbital selective CDW order- the metal remains with only short range
correlation despite its hybridization to the long range CDW layer. This
conclusion is supported by the data in Fig. 6, where sweeps of the scaled
$S^{\rm cdw}$ with $t_{3}$ at two values $\beta=20,\,10$ show a pair of
crossings. For $\beta=20$ these occur at $t_{3}\sim 1.3$ and $t_{3}\sim 3.2$
and for $\beta=10$ at $t_{3}\sim 1.33$ and $t_{3}\sim 3.1$. This analysis
indicates that long range CDW order exists in layer $l=-1$ only between these
values of $t_{3}$. In particular, for $t_{3}$ less than the lower critical
value, layer $l=-1$ is not in a CDW phase while layer $l=+1$ is. All crossings
are obtained with critical exponents, of the two-dimensional Ising model, as
expected since a $Z_{2}$ symmetry is being broken.

##  VI Spectral Functions and Double Occupancies

Having examined structure factors and inter-layer correlators, we now turn to
the spectral functions and the double occupancies, both of which provide
additional insight into the ground state properties. The layer-resolved
spectral functions $A_{l}(\omega)$, shown in Fig. 7, the many-body analog of
the single particle density of states, provide confirming evidence for the
Holstein interface phase diagram of Fig. 1(b). In the top row, for small
$t_{3}$, the Holstein layer $l=+1$ exhibits a CDW gap. The gap at $t_{3}=0$ is
large; hybridization with the metal produces peaks closer to $\omega=0$, but a
smaller gap remains. On the other hand the metal layer $l=-1$ has finite Fermi
surface spectral weight $A_{-1}(\omega=0)\neq 0$, thus showing the OSCDW. In
the middle row, for intermediate $t_{3}$, both layers have a gap, consistent
with the measurement of simultaneous long range CDW order. Finally, in the
bottom row, for large $t_{3}$, both layers have finite Fermi surface spectral
weight $A_{l}(\omega=0)\neq 0$ for $l=+1,-1$. The system is in the charge
singlet (charge liquid) phase.

We note that although the bilayer and interface geometries have many
properties in common at large $t_{3}$, their spectral functions are different.
There is a gap in the bilayer case, but not for the interface. We have
verified, with separate exact diagonalization calculations, that for dimers
(i.e. $t_{3}\gg t$) with $\lambda_{+1}=\lambda_{-1}$ one finds $A(\omega)$ is
gapped, while when $\lambda_{+1}$ is nonzero $\lambda_{-1}=0$ one reproduces
the behavior shown in Fig. 7(c,f).

Figure 7:  Spectral function $A(\omega)$ at $\beta=12$ for several $t_{3}$
values cutting across the Holstein-metal interface phase diagram of Fig. 1.
Top: Small $t_{3}$. Middle: Intermediate $t_{3}$. Bottom: Large $t_{3}$. Left
and right columns correspond to Holstein and metallic layers $l=+1$ and $l=-1$
respectively.

In a perfect CDW phase, half of the sites are doubly occupied and half are
empty, and ${\cal D}=0.5$. In the absence of interactions, $\lambda=0$, all
four site occupation possibilities
$|0\rangle,\,|\uparrow\rangle,\,|\downarrow\rangle,$ and
$|\uparrow\downarrow\rangle$, are equally likely and ${\cal D}=0.25$. This is
also the case in the charge singlet phase. Figure 8 shows ${\cal D}$ as a
function of $t_{3}$. Panel (a) is for the bilayer, where ${\cal D}_{+1}={\cal
D}_{-1}$, and panel (b) for the interface geometry where the two are
inequivalent. In both cases, ${\cal D}_{+1}$ is seen to evolve between the CDW
and singlet limits, although it never attains the value ${\cal D}=0.5$ owing
to the presence of quantum fluctuations. For the interface, ${\cal D}_{-1}$
begins at the uncorrelated value, ${\cal D}_{-1}=0.25$ at $t_{3}=0$ since
$\lambda_{-1}=0$. The double occupancy evolution is quite similar to that of
the structure factor, Fig. 2(a,c). However, since ${\cal D}$ is a local
observable, it exhibits less sharp features than $S^{\rm cdw}$ in the vicinity
of the QCP and thus only provides qualitative evidence for a cross-over
between those two phases. Besides that, a simple model which exhibits a layer-
dependent trivial CDW formation with only electronic degrees of freedom, a
bilayer ionic model, already displays this characteristic non-monotonicity
with growing hybridization (see Appendix C).

##  VII Discussion and conclusions

Figure 8:  Double occupancy ${\cal D}$ shown as a function of interlayer
hopping $t_{3}$ at low temperature $\beta=20$. ${\cal D}$ is a local
observable, and its value is the same for $L=8$ and $L=12$ to within the
symbol size; we show the latter. (a) Holstein bilayer; (b) Holstein-metal
interface. In (a) the two layers are equivalent and a single curve is shown.
In (b) the green curve shows data on layer $l=-1$, whereas the grey curve
represents layer $l=+1$.

We have generalized our existing understanding of the effect of
interlayer/interorbital hybridization $t_{3}$ on magnetic order driven by an
on-site electron-electron repulsion in the Hubbard model to charge order
originating in electron-phonon interactions in the Holstein model. The two
scenarios, although qualitatively related, are quite distinct in detail owing
to the lower symmetry of the CDW order parameter relative to the magnetic
case. Despite this difference, and its consequences such as the appearance of
charge order at finite temperature, the basic feature of the destruction of
long range order in the limit of large hybridization is shown still to occur.
Indeed, one remarkable conclusion of our work is that $t_{3}$ seems to restore
the degeneracy of pairing and charge correlations at the QCP.

Our most interesting observation is that the coexistence of CDW order on a
layer with non-zero electron-phonon coupling $\lambda$ with a metallic phase
on the $\lambda=0$ layer, which is trivially true at $t_{3}=0$, likely extends
out to finite $t_{3}$. This conclusion is based on the inability to scale the
charge correlations in the $\lambda=0$ layer unless $t_{3}\gtrsim 1.4$ (for
$\lambda_{+1}=0.2$).

The possibility that charge order takes place selectively parallels the known
occurrence of distinct Mott transitions in multi-orbital Hubbard models in the
coexistence of metallic and insulating behavior. The connection is, however,
not exact, since in principle a Mott transition might occur in the absence of
spontaneous symmetry breaking, whereas the insulating CDW phase here breaks
$Z_{2}$ symmetry. With that said, the Mott transition in its most common
incarnation, the square lattice Hubbard model, is always accompanied by long
range antiferromagnetic order. Thus our work does provide a close analog of
the case of orbital selective transitions in bands with differing electron-
electron interaction strengths.

###### Acknowledgements.

The work of Y.X.Z., C.H.F., and R.T.S. was supported by the grant DE‐SC0014671
funded by the U.S. Department of Energy, Office of Science. R.M. acknowledges
support from the NSFC Grants No. U1930402, No. 12050410263, No. 12111530010,
No. 11974039, and No. 12222401.

##  Appendix A Connection to Magnetic Language

In the repulsive 2D Hubbard model the dominant physics at half-filling on a
bipartite lattice is anti-ferromagnetic order, characterized by the operators,

$\displaystyle\hat{S}^{j}_{x}$
$\displaystyle=\frac{1}{2}\big{(}\,\hat{S}^{j}_{+}+\hat{S}^{j}_{-}\,\big{)}=\frac{1}{2}\big{(}\,\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}+\hat{c}^{\dagger}_{j\downarrow}\hat{c}^{\phantom{\dagger}}_{j\uparrow}\,\big{)}$
$\displaystyle\hat{S}^{j}_{y}$
$\displaystyle=\frac{1}{2i}\big{(}\,\hat{S}^{j}_{+}-\hat{S}^{j}_{-}\,\big{)}=\frac{1}{2i}\big{(}\,\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}-\hat{c}^{\dagger}_{j\downarrow}\hat{c}^{\phantom{\dagger}}_{j\uparrow}\,\big{)}$
$\displaystyle\hat{S}^{j}_{z}$
$\displaystyle=\frac{1}{2}\big{(}\,\hat{n}_{j\uparrow}-\hat{n}_{j\downarrow}\,\big{)}=\frac{1}{2}\big{(}\,\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\uparrow}-\hat{c}^{\dagger}_{j\downarrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}\,\big{)}$
(14)

From these relations, and as a consequence of the spin SU(2) symmetry of the
Hubbard model,

$\displaystyle\langle\hat{S}^{j}_{+}\hat{S}^{i}_{-}+\hat{S}^{j}_{-}\hat{S}^{i}_{+}\rangle$
$\displaystyle=4\,\langle\hat{S}^{j}_{z}\hat{S}^{i}_{z}\rangle$
$\displaystyle\langle\,\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}\hat{c}^{\dagger}_{i\downarrow}\hat{c}^{\phantom{\dagger}}_{i\uparrow}+\hat{c}^{\dagger}_{j\downarrow}\hat{c}^{\phantom{\dagger}}_{j\uparrow}\hat{c}^{\dagger}_{i\uparrow}\hat{c}^{\phantom{\dagger}}_{i\downarrow}\,\rangle$
$\displaystyle=$
$\displaystyle\langle\,\big{(}\,\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\uparrow}-\hat{c}^{\dagger}_{j\downarrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}\,\big{)}$
$\displaystyle\big{(}\,\hat{c}^{\dagger}_{i\uparrow}\hat{c}^{\phantom{\dagger}}_{i\uparrow}-\hat{c}^{\dagger}_{i\downarrow}\hat{c}^{\phantom{\dagger}}_{i\downarrow}\,\big{)}\,\rangle$
(15)

If we perform a particle-hole transformation to the down spin fermions,

$\displaystyle\hat{c}^{\dagger}_{j\downarrow}$
$\displaystyle\,\,\rightarrow\,\,(-1)^{j}\hat{c}^{\phantom{\dagger}}_{j\downarrow}$
$\displaystyle\hat{n}_{j\downarrow}\,\,$
$\displaystyle\rightarrow\,\,\big{(}\,1-\hat{n}_{j\downarrow}\,\big{)}$
$\displaystyle\hat{S}^{j}_{+}=\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\phantom{\dagger}}_{j\downarrow}$
$\displaystyle\,\,\rightarrow\,\,(-1)^{j}\hat{c}^{\dagger}_{j\uparrow}\hat{c}^{\dagger}_{j\downarrow}\equiv(-1)^{j}\hat{\Delta}^{\dagger}_{j}$
$\displaystyle\hat{S}^{j}_{z}=\frac{1}{2}\big{(}\,\hat{n}_{j\uparrow}-\hat{n}_{j\downarrow}\,\big{)}$
$\displaystyle\,\,\rightarrow\,\,\frac{1}{2}\big{(}\,\hat{n}_{j\uparrow}+\hat{n}_{j\downarrow}\,\big{)}\equiv\hat{n}_{j}$
(16)

we conclude that,

$\displaystyle-\,\langle\hat{\Delta}^{\dagger}_{j}\hat{\Delta}^{\phantom{\dagger}}_{i}+\hat{\Delta}^{\phantom{\dagger}}_{i}\hat{\Delta}^{\dagger}_{j}\rangle=\langle\big{(}\,\hat{n}_{j}-1)\,\big{)}\big{(}\,\hat{n}_{i}-1)\,\big{)}\rangle$
(17)

assuming that sites $i$ and $j$ are on opposite sublattices.

If, finally, assuming we are at half-filling, so that
$\langle\hat{n}_{j}\rangle=1$,

$\displaystyle-\,\langle\hat{\Delta}^{\dagger}_{j}\hat{\Delta}^{\phantom{\dagger}}_{i}+\hat{\Delta}^{\phantom{\dagger}}_{i}\hat{\Delta}^{\dagger}_{j}\rangle=\langle\hat{n}_{j}\hat{n}_{i}-1\rangle$
(18)

This shows that the two correlation functions of Eq. 5 are equal:
$p_{1,-1}=d_{1,-1}$. The merging of the two curves of Fig. 2(c,d) at a common
value reflects a restoration of an SU(2) symmetry of the Hubbard model. It is
interesting that this occurs even though the correlators are not in the
singlet limit of $-1/4$ (due to the fact that we are not in Holstein analog of
the large $U$ limit).

##  Appendix B Extracting the critical exponents

Figure 9: Contour plot of the sum of squared residuals of the least-squares
fit $S(\nu,\gamma)$ of the scaled data for the CDW structure factor $S^{\rm
cdw}$ of the Holstein bilayer at $t_{3}=2$ (see Fig. 3 for the original data).
A 16-th order polynomial is used to fit the data set, and the critical inverse
temperature used is $\beta_{c}=3.02$. The white marker denotes the minimum
$S(\nu,\gamma)$ which occurs at $\nu=0.95$ and $\gamma=1.7$, agreeing with the
2D Ising exponents.  Figure 10: Contour plot of the sum of squared residuals
of the least-squares fit $S(\nu,\gamma)$ of the scaled data for the CDW
structure factor $S^{\rm cdw}$ of the Holstein interface, top row at
$t_{3}=2$, bottom row at $t_{3}=1$, left column: layer $l=+1$, right column:
layer $l=-1$. (See Fig. 4 for the original data.) A 16-th order polynomial is
used to fit the data set, and the critical inverse temperature used are listed
in Fig. 5. The white marker denotes the minimum $S(\nu,\gamma)$ in the
displayed range of $\nu$ and $\gamma$. In panels (a,b,c), minimum residuals
are at $(\nu=0.95,\gamma=1.69)$, $(\nu=1.03,\gamma=1.76)$ and
$(\nu=0.89,\gamma=1.66)$ respectively, whereas in panel (d) the minimum occurs
at $\nu=0.83,\gamma=1.11$. Unlike panels (a,b,c), the $t_{3}=1$ metallic layer
of panel (d) does not fit the known Ising exponents, suggesting the absence of
order in this case.

We argue that owing to the symmetry of the order parameter of the CDW phase,
we make use of the critical exponents pertaining to the 2D Ising universality
class in order to simplify the FSS of the CDW structure factor $S^{\rm cdw}$.
Here, we justify this choice by quantitatively extracting the best set of
exponents $\nu$ and $\gamma$ that scales the curves according to the
functional form in Eq. 13. We start by using the scaled data in a large range
of critical exponents; subsequently, for each pair of $(\nu,\gamma)$, we
proceed with a high-order polynomial fitting of the scaled data, storing the
residual $S(\nu,\gamma)$ of the fitting procedure. The set of exponents that
minimizes $S(\nu,\gamma)$ is taken as those that characterize the transition.
The rationale is that if the dataset is well collapsed for a given
$(\nu,\gamma)$, a high-order polynomial fit (with the number of degrees of
freedom much smaller than the number of data points) will turn out to have a
fairly small error.

Using this procedure, we show in Fig. 9 the contour plot $S(\nu,\gamma)$ for
the data corresponding to the Holstein bilayer at $t_{3}=2$. The minimum
residual is indicated by the white marker. By observing its variation with
slightly different critical inverse temperature $\beta_{c}$, and different
polynomial orders used in the fits, we estimate $\nu=0.95\pm 0.07$ and
$\gamma=1.7\pm 0.1$, in agreement the 2D Ising exponents $\nu_{\rm 2D\
Ising}=1$ and $\gamma_{\rm 2D\ Ising}=7/4$. The formation of the CDW phase
establishes long range order by breaking the $Z_{2}$ sublattice symmetry,
which is the 2D Ising universality class.

We apply the same analysis to the interface model, as shown in Fig. 10 (See
Fig. 4 and Fig. 5 for original data). Both layers at $t_{3}=2$ as well as
layer $l=+1$ at $t_{3}=1$ minimize the residual close to the 2D Ising critical
exponents. However, in the case of layer $l=-1$ at $t_{3}=1$, the minimum
residue is located at $\nu=0.83,\gamma=1.11$ which is far from the 2D Ising
exponents. This supports our conclusion that there is no phase transition at
$t_{3}=1$ in layer $l=-1$ since a transition into a phase with long range CDW
order necessarily breaks the $Z_{2}$ symmetry and must be in the Ising
universality class.

##  Appendix C Induced CDW in Ionic Hubbard Model

We can get additional insight into the Holstein interface by considering the
following noninteracting, spinless, tight binding Hamiltonian,

$\displaystyle\mathcal{\hat{H}_{\rm BI-M}}=$ $\displaystyle-t\sum_{\langle
ij\rangle,l}\big{(}\,\hat{c}^{\dagger}_{i,l}\hat{c}^{\phantom{\dagger}}_{j,l}+{\rm
h.c.}\big{)}+\delta\sum_{i}(-1)^{i}\hat{n}_{i,1}$ $\displaystyle-
t_{3}\sum_{j}\big{(}\,\hat{c}^{\dagger}_{j,1}\hat{c}^{\phantom{\dagger}}_{j,-1}+{\rm
h.c.}\big{)}.$ (19)

Equation 19 describes two bands, labeled by $l=\pm 1$, each with hopping $t$
on a 2D square lattice, which are hybridized with each other by $t_{3}$. Band
$l=-1$ is metallic. At $t_{3}=0$ it has the usual 2D dispersion relation
$\epsilon(k)=-2t\,\big{(}\,{\rm cos}\,k_{x}+{\rm cos}\,k_{y}\,\big{)}$. Band
$l=+1$ is made insulating by the staggered potential $\delta$, so that at
$t_{3}=0$ its dispersion relation has two branches,
$E_{\pm}(k)=\pm\sqrt{\epsilon(k)^{2}+\delta^{2}}$. Both bands of Eq. 19 are
half-filled (the chemical potential $\mu=0$).

Figure 11:  Solution of the tight binding Hamiltonian, Eq. 19. (a) Solid
curves Occupations $\rho$ on the $+\delta$ sites of the insulating band as
functions of the magnitude of the staggered potential $\delta$; dashed curves:
Occupations on the partner sites in the metallic band to which those $+\delta$
sites are hybridized by $t_{3}$. (b) Solid Curve: Occupations on the $+\delta$
sites of the insulating band as a function of interlayer hybridization
$t_{3}$. Dashed curves: Occupations on the partner sites in the metallic band
to which those $+\delta$ sites are hybridized by $t_{3}$. The staggered
potential in this case is $\delta=1.25$. In both panels the linear lattice
size and the inverse temperature are $L=12$ and $\beta=4$, respectively.

In addition to inducing a band gap $2\delta$ in layer $l=+1$, the staggered
potential also creates a CDW phase, with low occupancy $+\delta$ and high
occupancy $-\delta$ sites. Here, the CDW order is trivial, in the sense of
being induced by an external field, as opposed to arising spontaneously in a
symmetric Hamiltonian like the Holstein model. Nevertheless we can still
examine how this ‘artificial’ CDW in layer $l=+1$ affects the site occupations
in the metallic band $l=-1$.

Figure 11(a) gives the occupations on the $+\delta$ sites of band $l=+1$ as
functions of $\delta$ for different $t_{3}$. As $\delta$ grows, the occupation
of the high energy sites in layer $l=+1$, which are directly coupled to the
staggered field, get increasingly small (solid curves). In turn, the
occupations of the partner sites on layer $l=-1$ which are not coupled to
$\delta$ are also shifted from $\rho=\frac{1}{2}$. These occupations increase
in order to take advantage of the decrease in the Pauli blocking. What is
interesting in the context of the simulations of the Holstein bilayers in the
main part of this paper is that, while the layer $l=+1$ occupations steadily
return to half-filling with increasing $t_{3}$, the evolution of the layer
$l=-1$ occupations is instead non-monotonic. The deviations of the occupations
from half-filling first grow with $t_{3}$, but then shrink.

This non-monotonicity is seen more clearly in Fig. 11(b) which plots similar
occupations as a function of $t_{3}$ for a fixed $\delta$. The maximum at
intermediate $t_{3}\sim 2.28$ is reminiscent of the behavior of Fig. 2(c),
which similarly shows a maximum in the induced CDW order at intermediate
$t_{3}$ in the metallic layer of the Holstein interface model. Indeed, the
agreement between the values of $t_{3}$ at which the induced order is maximal
is remarkably quantitative. To within error bars, the positions of the maxima
are the same, although the fall-off at large $t_{3}$ is more gradual in the
BI-Metal interface case.

## References

  * Bednorz and Müller [1988] J. G. Bednorz and K. A. Müller, Perovskite-type oxides—the new approach to high-${T}_{c}$ superconductivity, Rev. Mod. Phys. 60, 585 (1988).
  * Kastner  _et al._ [1998] M. Kastner, R. Birgeneau, G. Shirane, and Y. Endoh, Magnetic, transport, and optical properties of monolayer copper oxides, Rev. Mod. Phys. 70, 897 (1998).
  * Timusk and Statt [1999] T. Timusk and B. Statt, The pseudogap in high-temperature superconductors: an experimental survey, Rep. Prog. Phys. 62, 61 (1999).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * Vojta [2009] M. Vojta, Lattice symmetry breaking in cuprate superconductors: stripes, nematics, and superconductivity, Advances in Physics 58, 699 (2009).
  * Lopes dos Santos  _et al._ [2007] J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene bilayer with a twist: Electronic structure, Phys. Rev. Lett. 99, 256802 (2007).
  * Bistritzer and MacDonald [2010] R. Bistritzer and A. H. MacDonald, Transport between twisted graphene layers, Phys. Rev. B 81, 245412 (2010).
  * de Laissardiére  _et al._ [2010] G. T. de Laissardiére, D. Mayou, and L. Magaud, Localization of Dirac electrons in rotated graphene bilayers, Nano Lett. 10, 804 (2010).
  * Xue  _et al._ [2011] J. Xue, J. Sanchez-Yamagishi, D. Bulmash, P. Jacquod, A. Deshpande, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, and B. J. LeRoy, Scanning tunnelling microscopy and spectroscopy of ultra-flat graphene on hexagonal boron nitride, Nat. Mater. 10, 282 (2011).
  * Lang  _et al._ [2012] T. C. Lang, Z. Y. Meng, M. M. Scherer, S. Uebelacker, F. F. Assaad, A. Muramatsu, C. Honerkamp, and S. Wessel, Antiferromagnetism in the Hubbard model on the Bernal-stacked honeycomb bilayer, Phys. Rev. Lett. 109, 126402 (2012).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Rep. Prog. Phys 76, 056503 (2013).
  * 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).
  * 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).
  * Vekic  _et al._ [1995] M. Vekic, J. Cannon, D. Scalapino, R. Scalettar, and R. Sugar, Competition between antiferromagnetic order and spin liquid behavior in the two-dimensional periodic Anderson model at half-filling, Phys. Rev. Lett. 74, 2367 (1995).
  * 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).
  * 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).
  * Doniach [1977] S. Doniach, The Kondo lattice and weak antiferromagnetism, Physica B+C 91, 231 (1977).
  * Hess  _et al._ [1993] D. Hess, P. Riseborough, and J. Smith, in  _Encyclopedia of Applied Physics_, Vol. 7, edited by G. Trigg (VCH, New York, 1993) p. 435.
  * Grewe and Steglich [1991] N. Grewe and F. Steglich, in  _Handbook of Physics and Chemistry of Rare Earths_, Vol. 14, edited by J. K. A. Gschneidner and L. L. Eyring (Elsevier, Amsterdam, 1991) p. 343.
  * Liebsch [2004] A. Liebsch, Single Mott transition in the multiorbital Hubbard model, Phys. Rev. B 70, 165103 (2004).
  * Knecht  _et al._ [2005] C. Knecht, N. Blümer, and P. Van Dongen, Orbital-selective Mott transitions in the anisotropic two-band Hubbard model at finite temperatures, Phys. Rev. B 72, 081103 (2005).
  * Bouadim  _et al._ [2009] K. Bouadim, G. G. Batrouni, and R. T. Scalettar, Determinant quantum Monte Carlo study of the orbitally selective Mott transition, Phys. Rev. Lett. 102, 226402 (2009).
  * Koga  _et al._ [2004] A. Koga, N. Kawakami, T. M. Rice, and M. Sigrist, Orbital-selective Mott transitions in the degenerate Hubbard model, Phys. Rev. Lett. 92, 216402 (2004).
  * De Leo  _et al._ [2008] L. De Leo, M. Civelli, and G. Kotliar, ${T}=0$ heavy-fermion quantum critical point as an orbital-selective Mott transition, Phys. Rev. Lett. 101, 256404 (2008).
  * de’ Medici  _et al._ [2009] L. de’ Medici, S. R. Hassan, M. Capone, and X. Dai, Orbital-selective Mott transition out of band degeneracy lifting, Phys. Rev. Lett. 102, 126401 (2009).
  * Song  _et al._ [2015] Z.-Y. Song, H. Lee, and Y.-Z. Zhang, Possible origin of orbital selective Mott transitions in iron-based superconductors and Ca2-xSrxRuO4, New J. Phys. 17, 033034 (2015).
  * Song  _et al._ [2017] Z.-Y. Song, X.-C. Jiang, H.-Q. Lin, and Y.-Z. Zhang, Distinct nature of orbital-selective Mott phases dominated by low-energy local spin fluctuations, Phys. Rev. B 96, 235119 (2017).
  * Herbrych  _et al._ [2019] J. Herbrych, J. Heverhagen, N. Patel, G. Alvarez, M. Daghofer, A. Moreo, and E. Dagotto, Novel magnetic block states in low-dimensional iron-based superconductors, Phys. Rev. Lett. 123, 027203 (2019).
  * Herbrych  _et al._ [2020] J. Herbrych, G. Alvarez, A. Moreo, and E. Dagotto, Block orbital-selective Mott insulators: a spin excitation analysis, Phys. Rev. B 102, 115134 (2020).
  * Pandey  _et al._ [2020] B. Pandey, L.-F. Lin, R. Soni, N. Kaushal, J. Herbrych, G. Alvarez, and E. Dagotto, Prediction of exotic magnetic states in the alkali-metal quasi-one-dimensional iron selenide compound Na2FeSe2, Phys. Rev. B 102, 035149 (2020).
  * Środa  _et al._ [2021] M. Środa, E. Dagotto, and J. Herbrych, Quantum magnetism of iron-based ladders: Blocks, spirals, and spin flux, Phys. Rev. B 104, 045128 (2021).
  * Holstein [1959] T. Holstein, Studies of polaron motion: Part I. the molecular-crystal model, Ann. Phys. (N. Y.) 8, 325 (1959).
  * Scalettar  _et al._ [1989] R. T. Scalettar, N. E. Bickers, and D. J. Scalapino, Competition of pairing and Peierls–charge-density-wave correlations in a two-dimensional electron-phonon model, Phys. Rev. B 40, 197 (1989).
  * Marsiglio [1990] F. Marsiglio, Pairing and charge-density-wave correlations in the Holstein model at half-filling, Phys. Rev. B 42, 2416 (1990).
  * Weber and Hohenadler [2018] M. Weber and M. Hohenadler, Two-dimensional Holstein-Hubbard model: Critical temperature, Ising universality, and bipolaron liquid, Phys. Rev. B 98, 085405 (2018).
  * Costa  _et al._ [2018] N. C. Costa, T. Blommel, W.-T. Chiu, G. Batrouni, and R. Scalettar, Phonon dispersion and the competition between pairing and charge order, Phys. Rev. Lett. 120, 187003 (2018).
  * Hohenadler and Batrouni [2019] M. Hohenadler and G. Batrouni, Dominant charge density wave correlations in the Holstein model on the half-filled square lattice, Phys. Rev. B 100, 165114 (2019).
  * Li  _et al._ [2019] Z.-X. Li, M. L. Cohen, and D.-H. Lee, Enhancement of superconductivity by frustrating the charge order, Phys. Rev. B 100, 245105 (2019).
  * Zhang  _et al._ [2019] Y.-X. Zhang, W.-T. Chiu, N. Costa, G. Batrouni, and R. Scalettar, Charge order in the Holstein model on a honeycomb lattice, Phys. Rev. Lett. 122, 077602 (2019).
  * Chen  _et al._ [2019] C. Chen, X. Y. Xu, Z. Y. Meng, and M. Hohenadler, Charge-density-wave transitions of Dirac fermions coupled to phonons, Phys. Rev. Lett. 122, 077601 (2019).
  * Xiao  _et al._ [2021] B. Xiao, N. C. Costa, E. Khatami, G. G. Batrouni, and R. T. Scalettar, Charge density wave and superconductivity in the disordered Holstein model, Phys. Rev. B 103, L060501 (2021).
  * Cohen-Stead  _et al._ [2019] B. Cohen-Stead, N. C. Costa, E. Khatami, and R. T. Scalettar, Effect of strain on charge density wave order in the Holstein model, Phys. Rev. B 100, 045125 (2019).
  * Zhang  _et al._ [2020] Y.-X. Zhang, H.-M. Guo, and R. T. Scalettar, Charge density wave order on a $\pi$-flux square lattice, Phys. Rev. B 101, 205139 (2020).
  * Bradley  _et al._ [2021] O. Bradley, G. G. Batrouni, and R. T. Scalettar, Superconductivity and charge density wave order in the two-dimensional Holstein model, Phys. Rev. B 103, 235104 (2021).
  * [45] Studies of attractive Hubbard models, the anti-adiabatic limit of the Holstein model considered here, within the Dynamical Cluster Approximation are in progress. Philip Dee, Steven Johnston and Thomas Maier, private communication, work in progress.
  * Feng and Scalettar [2020] C. Feng and R. T. Scalettar, Interplay of flat electronic bands with Holstein phonons, Phys. Rev. B 102, 235152 (2020).
  * Batrouni and Scalettar [2019] G. G. Batrouni and R. T. Scalettar, Langevin simulations of a long-range electron-phonon model, Phys. Rev. B 99, 035114 (2019).
  * Loh  _et al._ [1990] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Quantum critical points and the sign problem, Science 375, 418 (2022).
  * Götz  _et al._ [2022] A. Götz, S. Beyl, M. Hohenadler, and F. Assaad, Valence-bond solid to antiferromagnet transition in the two-dimensional Su-Schrieffer-Heeger model by Langevin dynamics, Physical Review B 105, 085151 (2022).
  * Feng  _et al._ [2020] C. Feng, H. Guo, and R. T. Scalettar, Charge density waves on a half-filled decorated honeycomb lattice, Phys. Rev. B 101, 205103 (2020).