SciPost Phys. 15, 110 (2023)

Electron addition spectral functions of low-density polaron liquids
Alberto Nocera1,2⋆ and Mona Berciu1,2,3†
1 Stewart Blusson Quantum Matter Institute, University of British Columbia,
Vancouver, British Columbia, V6T 1Z4 Canada
2 Department of Physics Astronomy, University of British Columbia,
Vancouver, British Columbia, Canada V6T 1Z1
3 Leibniz Institute for Solid State and Materials Research (IFW) Dresden,
Helmholtzstrasse 20, 01069 Dresden, Germany⋆ alberto.nocera@ubc.ca ,† berciu@phas.ubc.ca

Abstract
Spectral functions are important quantities that contain a wealth of information about
the quasiparticles of a system, and that can also be measured experimentally. For systems with electron-phonon coupling, good approximations for the spectral function areavailable only in the Migdal limit (at Fermi energies much larger than the typical phonon
frequency, E F ≫ Ω, requiring a large carrier concentration x ) and in the single polaron
limit (at x = 0). Here we show that the region with x ≪ 1 (E F < Ω) can also be reliably investigated with the Momentum Average (MA) variational approximation, which
essentially describes the formation of a polaron above an inert Fermi sea. Specifically,
we show that for the one-dimensional spinless Holstein model, the MA spectral functions
compare favorably with those calculated using variationally exact density matrix renormalization group simulations (DMRG) evaluated directly in frequency-space, so long as
x < 0.1 and the adiabaticity ratio Ω/t > 0.5. Unlike in the Migdal limit, here ‘polaronic
physics’ emerges already at moderate couplings. The relevance of these results for a
spinful low-x metal is also discussed.
Copyright A. Nocera and M. Berciu.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 06-04-2023
Accepted 31-07-2023
Check forupdatesPublished 22-09-2023
doi:10.21468/SciPostPhys.15.3.110

Contents
1 Introduction

2

2 Model and results

3

3 Discussion
3.1 Methods
3.1.1 Momentum average approximation
3.1.2 Density matrix renormalization group

6
8
8
9

A Computational details to reproduce the DMRG results

10

References

13

1

SciPost Phys. 15, 110 (2023)

1

Introduction

A charged carrier polarizes its surroundings, dressing itself with a cloud of bosonic excitationsof the host crystal including phonons and charge, spin and/or orbital electronic fluctuations.
The properties of the resulting quasiparticle (when one is well defined) are of central practical interest because they play a key role in controlling the macroscopic behaviour of materials.
Solving such problems is also of tremendous fundamental interest because they provide a challenge for any formalism aiming to span weak, intermediate and strong coupling regimes—theirstudy has already led to many important developments in many-body physics and quantum
field theory.
An important subset of such problems is the study of polarons, which are the quasiparticles of systems with electron-phonon coupling. To calculate their spectral weights, which
can be directly measured with angle resolved photoemission spectroscopy, requires finding
the electron-addition propagator [1]†GSNe ,G(k, z) = GSNe ckσ Ĝ(z + EGS,Ne )ckσ(1)where Ĝ(z) = [z − H]−1 is the resolvent of the Hamiltonian H of interest, Ne = x N is thenumber of electrons set by the electron concentration x and the number of unit cells N → ∞,†H|GSNe 〉 = EGS,Ne |GSNe 〉 is the ground-state with Ne electrons, and ckσ
creates an electronwith momentum k and spin σ. We note that an electron-removal propagator can be definedsimilarly. [1] For the example considered below, it is obtained from the electron-addition propagator by replacing ω → −ω, x → 1 − x.
Despite considerable efforts and some significant successes in the ∼ 90 years since Landau
opened this field, [2] good general knowledge about the spectral functions of such systems isavailable only in two small regions of the parameter space: (i) single polarons, and (ii) the
Migdal limit.
The study of single polarons, at Ne = 0, is simpler because |GSNe =0 〉 ≡ |0〉 is the phonon vacuum and EGS,Ne =0 = 0. Moreover, the absence of other carriers render polaron-polaron interactions irrelevant. Single polarons have been studied for a variety of electron-phonon couplings,
the most famous being the Holstein, [3, 4] Fröhlich [5] and Peierls/SSH [6–9] models. Many
theoretical methods were developed to calculate the propagator of Eq. (1), including perturbation theory, [10] semiclassical approximations, [2,11] dynamical mean-field theory, [12,13]
path integral techniques, [14, 15] and variational methods [16] including the momentum average (MA) approximation. [17–20] These were supplemented by a wide variety of computational techniques [21] such as variational exact diagonalization [22, 23] and various typesof quantum Monte Carlo simulations, [24–28] plus the density matrix renormalization group(DMRG) in one dimension. [29, 30] By and large, there is now a good understanding of thegeneric properties of single polarons of model Hamiltonians at all couplings from very weak tovery strong, and many reliable tools to investigate new models. [31] More recently, significant
efforts have been devoted to combining these various tools with first principles calculations
for an accurate description of single polarons in actual materials. [32–34]
Generalizing these methods to study the spectral function of polarons in systems with finite carrier concentration x is difficult, first because the ground state |GS〉Ne now describes
a complicated liquid of polarons or bipolarons (i.e. two carriers bound through exchange ofphonons) [23, 35] and might exhibit superconductivity [36] or other orders. [28] Second,
while the additional electron again creates its own cloud, it also exchanges phonons with theother (bi)polarons, giving rise to effective electron-electron interactions that become considerable at stronger electron-phonon couplings. These strong effective interactions turn the system
into a strongly correlated one even if the bare electron-electron interactions are weak.

2

SciPost Phys. 15, 110 (2023)

It is possible to express the sum of the electron-addition and electron-removal propagators as an infinite series of diagrams which in principle are fully known [10] but in practice
are difficult to evaluate efficiently if the electron-phonon coupling is strong. Even for weakto-moderate couplings, an efficient semi-analytical summation is only possible in the Migdal
limit, for carrier concentrations x sufficiently large so that the Fermi energy E F is much larger
than the typical phonon energy Ω. Migdal showed that vertex corrections can be ignored
when E F ≫ Ω, greatly simplifying the summation of the remaining diagrams. [37] Even thislimit, however, cannot be treated accurately for strong electron-phonon couplings λ ∼ 1 where‘polaronic effects’ become considerable. Furthermore, it is well established that Migdal’s approximation fails when E F < Ω. [38]In this Letter we advance the study of electron-addition spectral functions to this unexplored area of the parameter space: low carrier concentrations x ≪ 1 so that E F < Ω, atall couplings. The results presented here are limited to spinless electrons to avoid the possibility of bipolaron formation and/or superconductivity. [39] As discussed below, we expect
these results to also be relevant for spinful electrons [40–45] if the bare electron-electron
repulsion is much stronger than the phonon-mediated attraction, preventing bipolaron formation/superconductivity.
Our results are obtained using a generalization of the MA variational approximation tosmall x concentrations, which describes the formation of a polaron when a fermion is added
above an inert Fermi sea. [46] We present MA results in one-dimension (1D) so that we canuse density matrix renormalization group (DMRG) –with significant enhancements obtainedin Ref. [47] and described in the Methods– as an unbiased computational method to gauge
them. This comparison validates the accuracy of the MA fermion addition spectral weights forx < 0.1 and Ω/t > 0.5 at all values of λ. It is important to emphasize that MA is trivial togeneralize to dimensions D > 1 and its accuracy improves with increasing D, [18] therefore
it allows the accurate study of polarons in any D for small x. The results presented here arefor the Holstein model, but DMRG (in 1D) and MA have been used successfully to study otherelectron-phonon couplings. [20] Thus, our work opens the way to efficiently study polaron
properties for a wide variety of models in a new region of the parameter space. Furthermore,valuable lessons are learnt from analyzing the processes that turn out to be most relevant
in MA, and may provide the key to expanding polaron studies to even wider regions of theparameter space.

2

Model and results

As advertised, we study the one-dimensional spinless Holstein model:

H = He + Hph + Ve−ph ,(2)Pwhere He = −t n (cn† cn+1 + H.c.) describes nearest-neighbor hopping on a chain with latticePpconstant a = 1, and cn† = k e−ikn ck† / N creates an electron at site n. Here N → ∞ isthe number of sites and k ∈Phonons are described by anP(−π, π] is the crystal momentum.
Einstein model: Hph = Ω n bn† bn (we set }h = 1) where bn† creates a a phonon at site n.
P
The Holstein electron-phonon coupling is Ve−ph = g n cn† cn (bn† + bn ). In standard fashion, wecharacterize the effective strength of electron-phonon coupling via the dimensionless coupling
λ = g 2 /(2Ωt).
The polaron dispersion is defined by the lowest-energy feature in the spectral weightA(k, ω) = − π1 ImG(k, ω + iη) where G(k, z) is the spinless analog of the propagator in Eq.
(1) and η → 0 is an artificial broadening. Before discussing MA and the approximations it
comprises, it is worth first seeing how it performs. Figures 1 and 2 show the electron-addition
3

SciPost Phys. 15, 110 (2023)

Figure 1: Electron-addition spectral functions A(k, ω) for Ω = t and η = 0.05.
From left to right, columns show results for x = 0, 0.05, 0.10, 0.15, respectively, whilerows show λ = 0.25, 0.50, 1.00 from top to bottom, respectively. Each contour plotshows the DMRG result for k < 0 and the MA result for k > 0. For x = 0 (left column)
the single polaron band extends across the entire Brillouin zone for all λ, moving tolower energies and becoming narrower as λ increases. In contrast, for weak coupling
λ = 0.25 and at finite x there is only a kink at E F + Ω. For moderate and large
λ = 0.5, 1, however, we see an emerging polaron band across at all momenta |k| > k F .
The polaron bandwidth renormalization decreases with increasing x, pointing to amore mobile quasiparticle. The agreement between MA and DMRG becomes poorer
with increasing x but is reasonable for x ≤ 0.1.
spectra when Ω = t and 0.5t, respectively, for concentrations x ∈ [0, 0.15] and couplings
λ ≤ 1. We note that our calculations are for a canonical ensemble so the energy ω is on an’absolute’ scale where ω = −2t marks the ground-state of a single bare fermion (x = 0, g = 0).
The more customary grandcanonical results are obtained by shifting the curves upwards by
the corresponding chemical potential µ, so that fermion addition spectral weight appears onlyfor ω ≥ 0. We also note that because MA approximates |GSNe 〉 with an inert Fermi sea, here
one only observes finite electron-addition spectral weight for |k| ≥ k F (the dashed verticals
lines mark k F ). DMRG confirms that this is accurate for the polaron band, however at higher
energies DMRG finds finite (albeit small) fermion addition spectral weight below k F .
The most striking observations are: (i) while a complete polaron band at all |k| ∈ [k F , π]
is observed for the single polaron (x = 0) as expected, it also appears for moderate-to-large
λ ≥ 0.5 for finite x. By contrast, for finite x and small λ = 0.25, we instead observe the kinkat E F + Ω predicted by perturbation theory, and also well established in the Migdal limit. The
emergence of the complete polaron band at finite x is a signature of the ‘polaronic effects’, andhere it is observed for a much weaker λ = 0.5 than the expected λ ∼ 1. We further discuss thisissue below. (ii) The polaron bandwidth narrowing (signalling a heavier polaron at stronger
coupling, all else being equal) decreases with increasing x, signalling a lighter polaron atlarger x (all else being equal). Results for the MA predictions of the polaron effective mass at

4

SciPost Phys. 15, 110 (2023)

Figure 2: Electron-addition spectral function A(k, ω) for Ω = 0.5t, with everythingelse as in Fig. 1. The evolution of the polaron band with x and λ mirrors the resultsfrom Fig. 1, but the quantitative agreement between DMRG and MA is worse than in
Fig. 1. For λ = 1, the polaron band is well separated from the higher energy features
and the spectrum is very unlike that of a free electron at all energies. See text for
more discussions.
various x, g values, and further discussion, are provided in Ref. [46]. The significant variationof the effective polaron mass m∗ with x even in this narrow range of x values, shows that
ignoring this dependence on doping in phenomenological models is very questionable. (iii)For λ = 0.25, the higher energy spectrum with ω > E F + Ω roughly follows the free particle
dispersion with additional broadening due to scattering on free phonons. In contrast, forstrong coupling λ = 1 the higher energy spectrum is split into a series of ’replicas’ that are
totally unlike the free electron spectrum. This change is also a clear signature of the ‘polaronic
effects’.
The existence or absence of a polaron state at larger k is rather hard to ascertain from Figs.
1,2 due to the broadening η = 0.05 used (this value was needed to get a better convergence
for the DMRG results, because the Correction Vector method [48–50] used here is known tobe numerically ill defined for η → 0 [51]). In Fig. 3, we plot the correspondings MA results
obtained at k = 0.8π, x = 0.1, Ω = 1 for several smaller η; these results are representative forall k near the Brillouin zone edge. They clearly show that while a polaron state is absent inthis region for λ = 0.25, it has already emerged for λ = 0.5; further increasing λ pushes thepolaron away from the continuum and gives it more quasiparticle weight. The appearance,in the low-density limit, of the polaron band at significantly smaller couplings than λ ∼ 1
expected in the Migdal limit raises questions about how accurate the latter estimate is. We
note that some of this difference could be due to dimensionality effects, something that canbe explored with MA.

5

SciPost Phys. 15, 110 (2023)

Figure 3: Emergence of the high-momentum polaron with increasing λ: Spectral
weight A(k = 0.8π, ω) for λ = 0.25, 0.5 and 1, respectively, and for different η
when x = 0.1, Ω = 1. For η = 0.05 (the value used in Figs. 1,2) the low energy
feature looks quite similar in all three cases. The lower η results demonstrate thatfor λ = 0.25, this feature is the lower edge of a continuum, i.e. here there is nopolaron state. For λ = 0.5, this feature resolves into a sharp peak lying below a
continuum, demonstrating that a polaron has formed. For λ = 1 the feature is asingle sharp peak, showing the polaron well below higher energy states.

3

Discussion

Overall, the agreement between DMRG and MA is qualitatively reasonable at most shown
values, becoming quantitatively accurate for x < 0.1 and improving for larger Ω/t ratios; allthis is achieved with a trivial MA computational cost (the results for any panel from Figs. 1
or 2 take under one second to generate on a regular workstation). In contrast, for the DMRG
calculations, a single frequency spectral function calculation can take up to 4 hours to complete
on a cluster node (we used an Intel Xeon Gold 6130 CPU node, hyperthreading over 8 cores);trivial parallelization is used to compute multiple frequencies simultaneously.
To understand the lessons that can be learned from these results, we now turn to analyzing the approximations made in MA. The first major approximation in the finite-x MA isto replace the unknown |GSNe 〉 by the mean-field approximation |m f 〉Ne = |F S〉 ⊗ |0̃〉, where
Q†|F S〉 =
|k|<k F ck |0〉 is the spinless Fermi sea with k F = πx, and |0̃〉 is a coherent phononstate describing the uniform lattice distortion induced by the uniform carrier distribution (seeMethods for more details). [46] We only allow one phonon cloud to be created when the extra
fermion is added, thus our variational method essentially analyzes the existence and the nature of a single polaron above an inert Fermi sea. As such, one can think of it as the analog of
Cooper’s calculation for the binding of two electrons above an inert Fermi sea, in the presenceof a weak attraction. [52] In that problem, all the electrons bind into pairs so one should use
the BCS wavefunction [36] instead of the inert Fermi sea; nevertheless, Cooper’s calculation
gives the correct binding energy because the BCS state and the Fermi sea only differ within avery narrow strip of width Ω centered at E F ≫ Ω. Similarly, in our problem all fermions get
dressed into polarons and the true ground state is some complicated polaron liquid. Because
here Ω > E F , one would naively expect this polaron liquid to be totally different from an inert
Fermi sea, nevertheless we find that the inert Fermi sea approximation works surprisingly wellfor small x.
The reason for this unexpected outcome can be qualitatively understood as follows: consider a polaron with a total momentum k T . Its energy is lowered due to level repulsion between
the configuration ck† |0〉 with energy ε(k T ) = −2t cos(k T ), and the one-phonon continuumT†comprising the configurations ck† bq=k|0〉 with energies ε(k) + Ω. For x = 0 all the oneT −kphonon configurations are available and the hybridization is strongest with the lowest-energy

6

SciPost Phys. 15, 110 (2023)

0.5

(a)

Ω=t
0.4

λ = 0.25

N − (k)λ = 0.5x=0
x=0.05x=0.10x=0.15

0.3

(b)

Ω=t

(c)

Ω=t
λ=1

0.2

0.1

0.0

0.5

(d)

Ω = 0.5t
0.4

λ = 0.25(e)Ω = 0.5t
λ = 0.5(f )Ω = 0.5t
λ=1

N − (k)

0.3

0.2

0.1

0.0
-π/2

-π/4

0
k

π/4

π/2

-π/4

0
k

π/4

π/2

-π/4

0
k

π/4

π/2

Figure 4:
Ground-state momentum distributions from DMRG:
N − (k) = 〈GSNe |ck† ck |GSNe 〉 vs momentum k at various densities x, couplings λfor Ω/t = 0.5 and 1. By contrast, the ’inert Fermi sea’ approximation built into thisMA implies N − (k) = Θ(k F − |k|), i.e. a step-like function set by the Fermi momentum
k F = xπ.
ones with k ≈ 0, q ≈ k T . At finite x, however, the small |k| fermion states are already occupied due to the formation of other polarons. Hybridization with the one-phonon continuumis therefore partially inhibited, explaining why polarons in the polaron liquid are less dressedthan single polarons of the same momentum and hence, why the polaron liquid wavefunction is not as different from the inert Fermi sea as one would naively expect for E F < Ω. [53]This argument also explains the decreased dressing with increasing x of the polaron withk > k F , [46, 53] shown by Figs. 1 and 2.
DMRG can test the validity of this ‘inert Fermi sea’ approximation, by calculating
N − (k) = 〈GSNe |ck† ck |GSNe 〉, shown in Fig. 4 for representative couplings and densities. For
any x > 0 these distributions are peaked at k = 0 and decrease on a scale set by k F = xπ,
however, they are far from the Heaviside function Θ(k F − |k|) of an inert Fermi sea used byMA.
The second major approximation in this finite-x MA is to ignore phonon renormalization by
not allowing particle-hole pairs to be excited out of the Fermi sea through phonon absorption.
Interestingly, this approximation is also used in the Migdal limit, raising the question whether
it is justified everywhere in the parameter space. In fact, we can include these processes
relatively easily in MA, however for x < 0.15 we found that they have hardly any effect on thepolaron band (differences are hard to see on the scale of Figs. 1 and 2). Higher energy features
are more strongly affected and the agreement with DMRG improves, however not enough, in
7

SciPost Phys. 15, 110 (2023)

our opinion, to justify the increased computational cost. We will report these incremental
improvements elsewhere.
The third approximation is the same as for the x = 0 MA, namely the variational constraints
on the phonon configurations included in the polaron cloud. The results in Figs. 1 and 2 arefrom a one-site MA(2) approximation, [19] which allows the phonon cloud to extend over one
site (with arbitrary number of phonons at this site which can be at any distance from theadditional fermion); in addition, up to two free phonons not bound to the cloud can be placed
anywhere else. We have also implemented a three-site MA(0) approximation, [20] where thephonon cloud is allowed to extend over three consecutive sites and there are no free phonons.
The results are qualitatively similar and the quantitative differences for the polaron bands are
not significant enough to justify the increased computational cost of the bigger variational
space.
Alltogether, this analysis reveals that the main step necessary to improve the agreement forlarger x is to find a better description of the |GSNe 〉 polaron liquid. In other words, we need tofind the polaron liquid analog to the BCS wavefunction. Nevertheless, we emphasize that thefinite-x MA presented here already works well for small x, is easy to implement and is very
efficient, and can be straightforwardly generalized to other types of bosons and of couplings
as well as to higher D, becoming more accurate with increasing D. It therefore allows us to
explore with confidence a part of the parameter space that was not studied before, and whichis relevant for many quantum materials that show interesting behavior at low doping, whether
it is in a nearly empty band (small x) or a nearly full band (small 1 − x).
Finally, we note that if the phonon-mediated attraction is small compared to the bareelectron-electron attraction so that bipolarons are unstable, this finite-x MA carries essentially without change to spinful fermions. Furthermore, because MA has been generalized to
study single bipolarons, [54] we can also study the effect of electron-phonon coupling on two
fermions above an inert Fermi sea, generalizing Cooper’s calculation away from weak couplingand Ω ≪ E F , and possibly opening the way to understand finite-x bipolaron liquids. This work
is under way.

3.1

Methods

3.1.1

Momentum average approximation

The finite-x MA generalization for both electron addition and electron removal propagators ispresented in Ref. [46], up to one important correction which we review here. This has to dowith the definition of the dressed phonon operators defining the new ’vacuum’ Bi |0̃〉 = 0 (forthe original phonons, |0̃〉 is a coherent state describing the mean-field uniform lattice distortion
induced by the uniform density of fermions). Instead of defining them as Bi = bi + g x/Ω,
[46] one needs to use Bi = bi + g N̂e /(ΩN ) where N̂e is the operator for the total number of
fermions. In the ground-state N̂e → Ne = N x and the two definitions agree, however whenthe extra fermion is added, N̂e → Ne + 1. This subtle difference is important when considering
the contribution from the mean-field shift −g 2 N̂e2 /(ΩN ) to the difference H − EGS,Ne in the
resolvent. In Ref. [46] these energy shifts cancelled out. The proper definition used here leads
to a finite difference in the thermodynamic limit −g 2 (Ne +1)2 /(ΩN )+ g 2 Ne2 /(ΩN ) = −2g 2 x/Ω
which has a significant effect on the location of the spectral weight for the larger x, λ values.
This correction is trivially accounted for by replacing z → z + 2g 2 x/Ω everywhere in the MA
expressions provided in Ref. [46].

8

SciPost Phys. 15, 110 (2023)

3.1.2

Density matrix renormalization group

We use the DMRG [55] as implemented in the DMRG++ software [56] at zero temperature,
to calculate the spectral functions A(k, ω) using the recently developed DMRG root-N Krylov
correction-vector approach [47].
In a nutshell, the root-N Krylov DMRG method evaluates the “Correction-Vector” [48,49] |CV x c 〉 = [ω − H + EGS,Ne + iη]−1 c †x |GSNe 〉 by applying N times the root-N progatorc[ω− H + EGS,Ne + iη]−1/N to the initial vector c †x |GSNe 〉, using at each step the standard DMRG
cCorrection Vector Krylov algorithm introduced in Ref. [50]. In this work, we use finite-size
chains with open boundary conditions (OBC) and employ the center-site approximation, in
which the additional fermion is added at the center site x c = N /2. This approximation reduces
the computational cost by order of N and becomes exact in the thermodynamic limit, although
it introduces “ringing” artifacts in the spatial Fourier Transform in small finite systems. Indeed,once the Correction Vector is computed, one extracts Oi,x c (ω) = − π1 Im[〈GSNe |ci |CV x c 〉] for all
PN
sites i. These are then used to calculate A(k, ω) = i=1 cos[k(ri − rc )]Oi,x c (ω). As seen in
figs. 1-2, in systems with N = 80 lattice sites artifacts from the use of OBC and the center site
approximation are minor.
In Ref. [47] it was shown that, at sufficiently large N , the DMRG root-N Krylov method
improves both the computational speed and frequency resolution (signal-to-noise ratio) of
the spectal functions compared to the standard Correction-Vector DMRG method [48–50]. In
particular, the root-N Krylov method is especially useful at large target frequencies with respect
to the Fermi energy, where large bond dimensions (or DMRG states) are typically needed.
These improvements were essential to obtain the DMRG results shown here. We refer thereader to Ref. [47] for more details.
In this work, we used the standard Fock basis representation of the phonon degrees of
freedom, utilizing up to 12 phonon states to represent the local phonon Hilbert space. We
therefore did not need to use the more sophisticated local phonon optimization methods [29,
57–63] or the recently developed projected purification DMRG method. [64, 65] Numericalresults were converged with respect to the bond dimension m. A maximum m = 1000 (and aminimum mmin = 24) provides convergence with a truncation error smaller than 10−6 for thefrequency dependent calculations. For the root-N Correction Vector Krylov calculations the
choice N = 20 has shown the best compromise in terms of moderately large bond dimension
required and computational speed. Finally, we set the the Krylov space tridiagonalization
error to εTridiag = 10−9 in order to avoid the proliferation of Krylov vectors (and thus Lanczos
iterations), and their reorthogonalizations. Sec. A provides computational details (input filefor the DMRG++ code) to reproduce the DMRG data shown in this work.

Data availability
The instructions to build the input scripts for the DMRG++ package to reproduce the DMRG
results can be found in Sec. A. The DMRG data/scripts to reproduce our figures is available asa public data set at [66]. Raw DMRG data files will be made available upon request.

Code availability
The numerical results reported in this work were obtained with DMRG++ versions 6.05 andPsimagLite versions 3.04. The DMRG++ computer program [56] is available at https://github.
com/g1257/dmrgpp.git, see Sec. A for more details.

9

SciPost Phys. 15, 110 (2023)

Acknowledgements
We thank A. Feiguin, S. Johnston, N. Prokof’ev, G. Sawatzky, S. Paeckel, and J. Sous for useful comments and suggestions. We acknowledge support from the Max Planck-UBC-UTokyo
Center for Quantum Materials and Canada First Research Excellence Fund (CFREF) Quantum
Materials and Future Technologies Program of the Stewart Blusson Quantum Matter Institute(SBQMI), and the Natural Sciences and Engineering Research Council of Canada (NSERC).
This work used computational resources and services provided by the Advanced Research Computing at the University of British Columbia. MB acknowledges the hospitality of the Leibniz
Institute for Solid State and Materials Research (IFW) Dresden, where part of this work was
carried out.

Author contributions
M. Berciu conceived the project and performed the MA calculations. A. Nocera performed the
DMRG calculations. M. Berciu and A. Nocera wrote the manuscript.

Competing interests
The authors declare no competing interests.
Correspondence and requests for materials should be addressed to A. Nocera and/or M.
Berciu.

A

Computational details to reproduce the DMRG results

Here we provide instructions on how to reproduce the DMRG results used in the main text.
The results reported in this work were obtained with DMRG++ versions 6.05 and PsimagLite
versions 3.04. The DMRG++ computer program [56] can be obtained with:

git clone https://github.com/g1257/dmrgpp.git
git clone https://github.com/g1257/PsimagLite.git
The main dependencies of the code are BOOST and HDF5 libraries. To compile the program:cd PsimagLite/lib; perl configure.pl; makecd ../../dmrgpp/src; perl configure.pl; make
The DMRG++ documentation can be found athttps://g1257.github.io/dmrgPlusPlus/manual.html or can be obtained by doing
cd dmrgpp/doc; make manual.pdf. In the description of the DMRG++ inputs below, we
follow very closely the description in the supplemental material of Ref. [47], where similarcalculations were performed.
The spectral function results for the 1D Holstein model for L = 80 sites, N = 8 electrons (thus x = 0.1), Ω = t and λ = 0.25 can be reproduced as follows. We first run./dmrg -f inputGS.ain -p 12 to obtain the ground state wave-function and ground stateenergy with 12 digit precision using the -p 12 option. The inputGS.ain has the form (thisis provided at [66])10SciPost Phys. 15, 110 (2023)##Ainur1.0
TotalNumberOfSites=160;
NumberOfTerms=3;### \sum_{<i,j>} -t * (c^\dag_i c_j + h.c.)
gt0:DegreesOfFreedom=1;
gt0:GeometryKind="ladder";gt0:GeometryOptions="ConstantValues";
gt0:dir0:Connectors=[1.0];
gt0:dir1:Connectors=[0.0];
gt0:LadderLeg=2;### bosonic hopping \sum_{<i,j>} -t_B * (b^\dag_i b_j + h.c.)gt1:DegreesOfFreedom=1;
gt1:GeometryKind="ladder";gt1:GeometryOptions="ConstantValues";
gt1:dir0:Connectors=[0.0];
gt1:dir1:Connectors=[0.0];
gt1:LadderLeg=2;### electron-phonon interaction \sum_i g*c^\dag_i c_i * (b^\dag_i+b_i)### with g=sqrt(2*Omega*lambda)
gt2:DegreesOfFreedom=1;
gt2:GeometryKind="ladder";gt2:GeometryOptions="ConstantValues";
gt2:dir0:Connectors=[0.0];
gt2:dir1:Connectors=[0.707];
gt2:LadderLeg=2;
Model="HolsteinSpinlessThin";
SolverOptions="twositedmrg,CorrectionTargeting,
vectorwithoffsets,useComplex";InfiniteLoopKeptStates=24;
FiniteLoops=[[79, 1000, 0],[-158, 1000, 0],[158, 1000, 0]];# Keep a maximum of 1000 states, but allow SVD truncation with# tolerance 1e-12 and minimum states equal to 24
TruncationTolerance="1e-12,24";# Symmetry sector for ground state N_e = 8
TargetElectronsTotal=8;# Associated with CorrectionTargeting: noise strength added to# reduced density matrix to avoid the algorith gets stuck;
CorrectionA=0.01;
OutputFile="dataGS_L80_N8_nph8";
The next step is to calculate dynamics for the A(k, ω) spectral function using the saved groundstate as an input. It is convenient to do the dynamics run in a subdirectory Aqw, so create
this directory first, and then add/modify the following lines in inputAqw.ado (this input isprovided at [66])# The finite loops now start from the final loop of the gs calculation.
# Total number of finite loops equal to N+2, here N=6
FiniteLoops=[[-158, 1000, 2],[158, 1000, 2],[-158, 1000, 2],[158, 1000, 2],[-158, 1000, 2],[158, 1000, 2],[-158, 1000, 2],[158, 1000, 2]];11

SciPost Phys. 15, 110 (2023)# Keep a maximum of 1000 states, but allow SVD truncation with# tolerance 1e-6 and minimum states equal to 24
TruncationTolerance="1e-6,24";# The exponent in the root-N CV method
CVnForFraction=6;# Tolerance for Tridiagonal Decomposition of the effective Hamiltonian
TridiagonalEps=1e-12;# Solver options should appear on one line, here we have two lines because# of formatting purposes
SolverOptions="useComplex,twositedmrg,vectorwithoffsets,
TargetingCVEvolution,restart,fixLegacyBugs,minimizeDisk";# RestartFilename is the name of the GS .hd5 file (extension is not needed)
RestartFilename="../dataGS_L80_N8_nph8";# The weight of the g.s. in the density matrix
GsWeight=0.1;# Legacy, set to 0
CorrectionA=0;# Fermion spectra has sign changes in denominator.
# For boson operators (as in here) set it to 0
DynamicDmrgType=0;# The site(s) where to apply the operator below. Here it is the center site.
TSPSites=[79];# The delay in loop units before applying the operator. Set to 0
TSPLoops=[0];# If more than one operator is to be applied, how they should be combined.
# Irrelevant if only one operator is applied, as is the case here.
TSPProductOrSum="sum";# Sets the number of sweeps to 1 before advancing in "time"=1/N
TSPAdvanceEach=78;# How the operator to be applied will be specified
string TSPOp0:TSPOperator=expression;# The operator expression to apply the c^\dagger operator on the center site
string TSPOp0:OperatorExpression="c?0’";# How is the freq. given in the denominator (Matsubara is the other option)
CorrectionVectorFreqType="Real";# This is a dollarized input, so the# omega will change from input to input.
CorrectionVectorOmega=$omega;# The broadening for the spectrum in omega + i*eta
CorrectionVectorEta=0.05;# The algorithm
CorrectionVectorAlgorithm="Krylov";#The labels below are ONLY read by manyOmegas.pl script# How many inputs files to create#OmegaTotal=40# Which one is the first omega value#OmegaBegin=-4.0# Which is the "step" in omega#OmegaStep=0.1# Because the script will also be creating the batches,# indicate what to measure in the batches
12

SciPost Phys. 15, 110 (2023)#Observable=c?0
Notice that the main change with respect of a standard CV method input is given by the option TargetingCVEvolution in the SolverOptions instead of CorrectionVectorTargeting,and the addition of the line CVnForFraction=6; We note also that the number of finite loops
must be at least equal to the number equal to the exponent in the root-N CV method. As inthe standard CV approach, all individual inputs (one per ω in the correction vector approach)can be generated and submitted using the manyOmegas.pl script which can be found in the
dmrgpp/src/script folder (but also provided at [66]):

perl manyOmegas.pl inputAqw.ado batchTemplate.pbs <test/submit>.
It is recommended to run with test first to verify correctness, before running with submit.
Depending on the machine and scheduler, the BatchTemplate can be e.g. a PBS or SLURM
script. The key is that it contains a line ./dmrg -f $$input "<gs|$$obs|P1>" -p 12
which allows manyOmegas.pl to fill in the appropriate input for each generated job batch.
After all outputs have been generated,

perl myprocAkw.pl inputAqw.ado
can be used to process (this script is also provided at [66]) and generate a data file
outSpectrum.c?0.gnuplot ready to be plotted (using the Gnuplot software, for example).

References[1] A. Damascelli, Z. Hussain and Z.-X. Shen, Angle-resolved
studies of the cuprate superconductors, Rev. Mod. Phys. 75,
doi:10.1103/RevModPhys.75.473.

photoemission
473 (2003),[2] L. D. Landau, Electron motion in crystal lattices, Phys. Z. Sowjetunion 3, 664 (1933),
doi:10.1016/b978-0-08-010586-4.50015-8.
[3] T. Holstein, Studies of polaron motion: Part I. The molecular-crystal model, Ann. Phys. 8,325 (1959), doi:10.1016/0003-4916(59)90002-8.
[4] T. Holstein, Studies of polaron motion: Part II. The “small” polaron, Ann. Phys. 8, 343(1959), doi:10.1016/0003-4916(59)90003-X.
[5] H. Fröhlich,
Electrons in latticedoi:10.1080/00018735400101213.

fields,

Adv.

Phys.

3,

325

(1954),[6] R. Peierls, Quantum theory of solids, Oxford University Press, Oxford, UK, ISBN
9780192670175 (1955), doi:10.1093/acprof:oso/9780198507819.001.0001.
[7] S. Barišić, J. Labbé and J. Friedel, Tight binding and transition-metal superconductivity,
Phys. Rev. Lett. 25, 919 (1970), doi:10.1103/PhysRevLett.25.919.
[8] W. P. Su, J. R. Schrieffer and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42,
1698 (1979), doi:10.1103/PhysRevLett.42.1698.
[9] A. J. Heeger, S. Kivelson, J. R. Schrieffer and W. -P. Su, Solitons in conducting polymers,
Rev. Mod. Phys. 60, 781 (1988), doi:10.1103/RevModPhys.60.781.
[10] G. D. Mahan, Many-particle physics, Springer, New York, USA, ISBN 9781475757149(2013), doi:10.1007/978-1-4757-5714-9.
13

SciPost Phys. 15, 110 (2023)[11] L. D. Landau and S. Pekar, Effective mass of a polaron, Zh. Eksp. Teor. Fiz. 18, 419 (1948),
doi:10.1016/b978-0-08-010586-4.50072-9.
[12] S. Ciuchi, F. de Pasquale, S. Fratini and D. Feinberg, Dynamical mean-field theory of thesmall polaron, Phys. Rev. B 56, 4494 (1997), doi:10.1103/PhysRevB.56.4494.
[13] P. Mitrić, V. Janković, N. Vukmirović and D. Tanasković, Spectral functions of the Holstein polaron: Exact and approximate solutions, Phys. Rev. Lett. 129, 096401 (2022),
doi:10.1103/PhysRevLett.129.096401.
[14] R. P. Feynman, Slow electrons in a polar crystal, Phys. Rev. 97, 660 (1955),
doi:10.1103/PhysRev.97.660.
[15] R. P. Feynman, R. W. Hellwarth, C. K. Iddings and P. M. Platzman, Mobility of slow electrons
in a polar crystal, Phys. Rev. 127, 1004 (1962), doi:10.1103/PhysRev.127.1004.
[16] Y. Toyozawa, Self-trapping of an electron by the acoustical mode of lattice vibration. I, Prog.
Theor. Phys. 26, 29 (1961), doi:10.1143/PTP.26.29.
[17] M. Berciu, Green’s function of a dressed particle, Phys. Rev. Lett. 97, 036402 (2006),
doi:10.1103/PhysRevLett.97.036402.
[18] G. L. Goodvin, M. Berciu and G. A. Sawatzky, Green’s function of the Holstein polaron,
Phys. Rev. B 74, 245104 (2006), doi:10.1103/PhysRevB.74.245104.
[19] M. Berciu and G. L. Goodvin, Systematic improvement of the momentum average approximation for the Green’s function of a Holstein polaron, Phys. Rev. B 76, 165109 (2007),
doi:10.1103/PhysRevB.76.165109.
[20] D. J. J. Marchand, G. De Filippis, V. Cataudella, M. Berciu, N. Nagaosa, N. V.
Prokof’ev, A. S. Mishchenko and P. C. E. Stamp, Sharp transition for single polarons in
the one-dimensional Su-Schrieffer-Heeger model, Phys. Rev. Lett. 105, 266605 (2010),
doi:10.1103/PhysRevLett.105.266605.
[21] H. Fehske and S. A. Trugman, Numerical solution of the Holstein polaron problem, in
Polarons in advanced materials, Springer, Dordrecht, Netherlands, ISBN 9781402063473(2007), doi:10.1007/978-1-4020-6348-0_10.
[22] J. Bonča, S. A. Trugman and I. Batistić, Holstein polaron, Phys. Rev. B 60, 1633 (1999),
doi:10.1103/PhysRevB.60.1633.
[23] J. Bonča and S. A. Trugman, Bipolarons in the extended Holstein Hubbard model, Phys.
Rev. B 64, 094507 (2001), doi:10.1103/PhysRevB.64.094507.
[24] N. V. Prokof’ev and B. V. Svistunov, Polaron problem by diagrammatic quantum Monte
Carlo, Phys. Rev. Lett. 81, 2514 (1998), doi:10.1103/PhysRevLett.81.2514.
[25] A. Mishchenko, N. Prokof’ev, A. Sakamoto and B. Svistunov, Diagrammatic quantum Monte Carlo study of the Fröhlich polaron, Phys. Rev. B 62, 6317 (2000),
doi:10.1103/PhysRevB.62.6317.
[26] 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),
doi:10.1103/PhysRevB.103.235104.

14

SciPost Phys. 15, 110 (2023)[27] B. Nosarzewski, E. W. Huang, P. M. Dee, I. Esterlis, B. Moritz, S. A. Kivelson, S. Johnston
and T. P. Devereaux, Superconductivity, charge density waves, and bipolarons in the Holsteinmodel, Phys. Rev. B 103, 235156 (2021), doi:10.1103/PhysRevB.103.235156.
[28] 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), doi:10.1103/PhysRevB.40.197.
[29] E. Jeckelmann and S. R. White, Density-matrix renormalization-group study ofthe polaron problem in the Holstein model, Phys. Rev. B 57, 6376 (1998),
doi:10.1103/PhysRevB.57.6376.
[30] C. Zhang, E. Jeckelmann and S. R. White, Dynamical properties of the one-dimensional
Holstein model, Phys. Rev. B 60, 14092 (1999), doi:10.1103/PhysRevB.60.14092.
[31] J. T. Devreese and A. S. Alexandrov, Fröhlich polaron and bipolaron: Recent developments,
Rep. Prog. Phys. 72, 066501 (2009), doi:10.1088/0034-4885/72/6/066501.
[32] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. 89, 015003(2017), doi:10.1103/RevModPhys.89.015003.
[33] W. H. Sio, C. Verdi, S. Poncé and F. Giustino, Polarons from first principles, without supercells, Phys. Rev. Lett. 122, 246403 (2019), doi:10.1103/PhysRevLett.122.246403.
[34] J. Lafuente-Bartolome, C. Lian, W. H. Sio, I. G. Gurtubay, A. Eiguren and F. Giustino,
Unified approach to polarons and phonon-induced band structure renormalization, Phys.
Rev. Lett. 129, 076402 (2022), doi:10.1103/PhysRevLett.129.076402.
[35] A. Nocera, J. Sous, A. E. Feiguin and M. Berciu, Bipolaron liquids atstrong Peierls electron-phonon couplings, Phys. Rev. B 104, L201109 (2021),
doi:10.1103/PhysRevB.104.L201109.
[36] J. Bardeen, L. N. Cooper and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. 108,
1175 (1957), doi:10.1103/PhysRev.108.1175.
[37] V. Galitskii and A. Migdal, Application of quantum field theory methods to the many bodyproblem, Sov. Phys. J. Exp. Theor. Phys. 7, 96 (1958).
[38] A. S. Mishchenko, I. S. Tupitsyn, N. Nagaosa and N. Prokof’ev, Fermi blockade of thestrong electron–phonon interaction in modelled optimally doped high temperature superconductors, Sci. Rep. 11, 9699 (2021), doi:10.1038/s41598-021-89059-w.
[39] E. Langmann, C. Hainzl, R. Seiringer and A. V. Balatsky, Obstructions to odd-frequency
superconductivity in Eliashberg theory, (arXiv preprint) doi:10.48550/arXiv.2207.01825.
[40] H. Matsueda, T. Tohyama and S. Maekawa, Electron-phonon coupling and spin-charge
separation in one-dimensional Mott insulators, Phys. Rev. B 74, 241103 (2006),
doi:10.1103/PhysRevB.74.241103.
[41] W.-Q. Ning, H. Zhao, C.-Q. Wu and H.-Q. Lin, Phonon effects on spin-charge separation in
one dimension, Phys. Rev. Lett. 96, 156402 (2006), doi:10.1103/PhysRevLett.96.156402.
[42] A. Nocera, M. Soltanieh-ha, C. A. Perroni, V. Cataudella and A. E. Feiguin, Interplay of charge, spin, and lattice degrees of freedom in the spectral propertiesof the one-dimensional Hubbard-Holstein model, Phys. Rev. B 90, 195134 (2014),
doi:10.1103/PhysRevB.90.195134.
15

SciPost Phys. 15, 110 (2023)[43] C. B. Mendl, E. A. Nowadnick, E. W. Huang, S. Johnston, B. Moritz and T. P. Devereaux,
Doping dependence of ordered phases and emergent quasiparticles in the doped HubbardHolstein model, Phys. Rev. B 96, 205141 (2017), doi:10.1103/PhysRevB.96.205141.
[44] E. A. Nowadnick, S. Johnston, B. Moritz and T. P. Devereaux, Renormalization of spectraby phase competition in the half-filled Hubbard-Holstein model, Phys. Rev. B 91, 165127(2015), doi:10.1103/PhysRevB.91.165127.
[45] Y. Wang, I. Esterlis, T. Shi, J. I. Cirac and E. Demler, Zero-temperature phases of the twodimensional Hubbard-Holstein model: A non-Gaussian exact diagonalization study, Phys.
Rev. Res.2, 043258 (2020), doi:10.1103/PhysRevResearch.2.043258.
[46] M. Berciu, Polarons in spinless metals — A variational solution, J. Phys. Mater. 5, 044002(2022), doi:10.1088/2515-7639/ac9086.
[47] A. Nocera and G. Alvarez, Root-n Krylov-space correction vectors for spectral functions
with the density matrix renormalization group, Phys. Rev. B 106, 205106 (2022),
doi:10.1103/PhysRevB.106.205106.
[48] T. D. Kühner and S. R. White, Dynamical correlation functions using the density matrix
renormalization group, Phys. Rev. B 60, 335 (1999), doi:10.1103/PhysRevB.60.335.
[49] E. Jeckelmann, Dynamical density-matrix renormalization-group method, Phys. Rev. B 66,
045114 (2002), doi:10.1103/PhysRevB.66.045114.
[50] A. Nocera and G. Alvarez, Spectral functions with the density matrix renormalizationgroup: Krylov-space approach for correction vectors, Phys. Rev. E 94, 053308 (2016),
doi:10.1103/PhysRevE.94.053308.
[51] M. Paech and E. Jeckelmann, Blind deconvolution of density-matrix renormalization-group
spectra, Phys. Rev. B 89, 195101 (2014), doi:10.1103/PhysRevB.89.195101.
[52] L. N. Cooper, Bound electron pairs in a degenerate Fermi gas, Phys. Rev. 104, 1189 (1956),
doi:10.1103/PhysRev.104.1189.
[53] A. S. Mishchenko, N. Nagaosa and N. Prokof’ev, Diagrammatic Monte Carlo
method for many-polaron problems, Phys. Rev. Lett. 113, 166402 (2014),
doi:10.1103/PhysRevLett.113.166402.
[54] J. Sous, M. Chakraborty, R. V. Krems and M. Berciu, Light bipolarons stabilized by Peierls electron-phonon coupling, Phys. Rev. Lett. 121, 247001 (2018),
doi:10.1103/PhysRevLett.121.247001.
[55] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev.
Lett. 69, 2863 (1992), doi:10.1103/PhysRevLett.69.2863.
[56] G. Alvarez, The density matrix renormalization group for strongly correlated electron systems: A generic implementation, Comp. Phys. Commun. 180, 1572 (2009),
doi:10.1016/j.cpc.2009.02.016.
[57] C. Zhang, E. Jeckelmann and S. R. White, Density matrix approach to local Hilbert space
reduction, Phys. Rev. Lett. 80, 2661 (1998), doi:10.1103/PhysRevLett.80.2661.
[58] C. Guo, A. Weichselbaum, J. von Delft and M. Vojta, Critical and strong-couplingphases in one- and two-bath spin-boson models, Phys. Rev. Lett. 108, 160401 (2012),
doi:10.1103/PhysRevLett.108.160401.
16

SciPost Phys. 15, 110 (2023)[59] C. Brockt, F. Dorfner, L. Vidmar, F. Heidrich-Meisner and E. Jeckelmann, Matrix-productstate method with a dynamical local basis optimization for bosonic systems out of equilibrium, Phys. Rev. B 92, 241106 (2015), doi:10.1103/PhysRevB.92.241106.
[60] D. Jansen, C. Jooss and F. Heidrich-Meisner, Charge density wave breakdown ina heterostructure with electron-phonon coupling, Phys. Rev. B 104, 195116 (2021),
doi:10.1103/PhysRevB.104.195116.
[61] J. Stolpp, T. Köhler, S. R. Manmana, E. Jeckelmann, F. Heidrich-Meisner and S. Paeckel,
Comparative study of state-of-the-art matrix-product-state methods for lattice models withlarge local Hilbert spaces without U(1) symmetry, Comput. Phys. Commun. 269, 108106(2021), doi:10.1016/j.cpc.2021.108106.
[62] J. Stolpp, J. Herbrych, F. Dorfner, E. Dagotto and F. Heidrich-Meisner, Charge-densitywave melting in the one-dimensional Holstein model, Phys. Rev. B 101, 035134 (2020),
doi:10.1103/PhysRevB.101.035134.
[63] D. Jansen, J. Bonča and F. Heidrich-Meisner, Finite-temperature optical conductivity with density-matrix renormalization group methods for the Holstein polaron and bipolaron with dispersive phonons, Phys. Rev. B 106, 155129 (2022),
doi:10.1103/PhysRevB.106.155129.
[64] T. Köhler, J. Stolpp and S. Paeckel, Efficient and flexible approach to simulate lowdimensional quantum lattice models with large local Hilbert spaces, SciPost Phys. 10, 058(2021), doi:10.21468/SciPostPhys.10.3.058.
[65] S. Mardazad, Y. Xu, X. Yang, M. Grundner, U. Schollwöck, H. Ma and S. Paeckel, Quantum
dynamics simulation of intramolecular singlet fission in covalently linked tetracene dimer,
J. Chem. Phys. 155, 194101 (2021), doi:10.1063/5.0068292.
[66] A. Nocera and M. Berciu, Electron addition spectral functions of low-density polaron liquids,
Zenodo (2023), doi:10.5281/zenodo.7790340.

17