# Quantum Adiabatic Doping for Atomic Fermi-Hubbard Quantum Simulations

Jue Nan  Hefei National Laboratory for Physical Sciences at Microscale and
Department of Modern Physics, University of Science and Technology of China,
Hefei, Anhui 230026, China.  State Key Laboratory of Surface Physics,
Institute of Nanoelectronics and Quantum Computing, and Department of Physics,
Fudan University, Shanghai 200433, China  Jian Lin  State Key Laboratory of
Surface Physics, Institute of Nanoelectronics and Quantum Computing, and
Department of Physics, Fudan University, Shanghai 200433, China  Yuchen Luo
State Key Laboratory of Surface Physics, Institute of Nanoelectronics and
Quantum Computing, and Department of Physics, Fudan University, Shanghai
200433, China  Bo Zhao<EMAIL_ADDRESS>Hefei National Laboratory for
Physical Sciences at Microscale and Department of Modern Physics, University
of Science and Technology of China, Hefei, Anhui 230026, China.  Shanghai
Branch,CAS Center for Excellence and Synergetic Innovation Center in Quantum
Information and Quantum Physics, University of Science and Technology of
China, Shanghai 201315, China.  Xiaopeng Li<EMAIL_ADDRESS>State
Key Laboratory of Surface Physics, Institute of Nanoelectronics and Quantum
Computing, and Department of Physics, Fudan University, Shanghai 200433, China
Shanghai Qi Zhi Institute, Xuhui District, Shanghai, 200032, China.

###### Abstract

There have been considerable research efforts devoted to quantum simulations
of Fermi-Hubbard model with ultracold atoms loaded in optical lattices. In
such experiments, the antiferromagnetically ordered quantum state has been
achieved at half filling in recent years. The atomic lattice away from half
filling is expected to host d-wave superconductivity, but its low temperature
phases have not been reached. In a recent work Lin  _et al._ (2019), we
proposed an approach of incommensurate quantum adiabatic doping, using quantum
adiabatic evolution of an incommensurate lattice for preparation of the highly
correlated many-body ground state of the doped Fermi-Hubbard model starting
from a unit-filling band insulator. Its feasibility has been demonstrated with
numerical simulations of the adiabatic preparation for certain incommensurate
particle-doping fractions, where the major problem to circumvent is the atomic
localization in the incommensurate lattice. Here we carry out a systematic
study of the quantum adiabatic doping for a wide range of doping fractions
from particle-doping to hole-doping, including both commensurate and
incommensurate cases. We find that there is still a localization-like slowing-
down problem at commensurate fillings, and that it becomes less harmful in the
hole-doped regime. With interactions, the adiabatic preparation is found to be
more efficient for that interaction effect destabilizes localization. For both
free and interacting cases, we find the adiabatic doping has better
performance in the hole-doped regime than the particle-doped regime. We also
study adiabatic doping starting from the half-filling Mott insulator, which is
found to be more efficient for certain filling fractions.

##  I Introduction

Ultracold atoms in optical lattices provide a fascinating platform for quantum
simulations of correlated many-body physics Bloch  _et al._ (2008); Li and Liu
(2016); Lewenstein  _et al._ (2007); Bloch (2018). Since the atomic tunneling
and interactions are both controllable in these systems, they have widely been
used to study quantum many-body phases and quantum phase transitions. One main
theme of quantum simulation with optical lattices is to investigate the low-
temperature phase diagram of the Fermi-Hubbard model Jaksch  _et al._ (1998);
Hofstetter  _et al._ (2002); Dutta  _et al._ (2015), and help uncover the
fundamental mechanism of high-temperature superconductivity Lee  _et al._
(2006).

Whether and how the d-wave superconductivity arises in the doped region
including both hole- and particle doped cases, in the repulsive Fermi-Hubbard
model has been attracting continuous research efforts Lee  _et al._ (2006);
Metzner  _et al._ (2012); Jiang and Devereaux (2019); Qin  _et al._ (2020),
but this remains an open question with no consensus reached Jiang and
Devereaux (2019); Qin  _et al._ (2020), one reason being that the numerical
simulations on classical computers meet fundamental challenges of
exponentially growing Hilbert space of the quantum many-body system. This
makes quantum simulations of the doped Fermi-Hubbard model very much demanded
as it is has potential to address the important question of existence of
d-wave superconductivity in the model. With the development of ultracold atom
experiments, the low-temperature antiferromagnetic phase has now been reached
at half filling Greif  _et al._ (2013); Hart  _et al._ (2015); Mazurenko  _et
al._ (2017). The doping of antiferromagnet with hole has been realized by
reducing the density of the trapped gases Mazurenko  _et al._ (2017). However,
it is difficult to cool down the system to a sufficiently low temperature to
enter the d-wave superconducting phase. How to optimally perform doping for
the atomic Fermi-Hubbard optical lattice system while maintaining low entropy
of the system demands more theoretical study.

A plausible method to maintain the quantum simulator at low entropy is to
perform the quantum adiabatic doping. The adiabatic quantum state preparation
has been widely applied in the quantum state engineering and quantum
simulation Trebst  _et al._ (2006); Sørensen  _et al._ (2010); Lubasch  _et
al._ (2011); Zhang and Duan (2013); Chiu  _et al._ (2018); Lin  _et al._
(2019); Sun  _et al._ (2020). For the fermions confined in the optical
lattice, a different filling factor is achieved by adiabatically converting
the lattice with one spatial period to a lattice with a different period.
Going through an adiabatic evolution from insulating states to the doped
regime, the system remains at the ground state of the instantaneous
Hamiltonian. This protocol has been studied to prepare Fermi-Hubbard
antiferromagnet insulating stateLubasch  _et al._ (2011); Chiu  _et al._
(2018), and also doped ground state Lin  _et al._ (2019) with incommensurate
fillings. For the incommensurate case, it has been found that the major
difficulty in carrying out the quantum adiabatic doping is from fermion
localization Lin  _et al._ (2019). The physics of fermion localization
occurring in the intermediate dynamics prevents efficient state preparation,
causing a problem of localization slowing down. As a generic solution of the
localization slowing down Basko  _et al._ (2006); Oganesyan and Huse (2007);
Pal and Huse (2010); Kondov  _et al._ (2015), the atomic interaction is
introduced and found to improve the preparation efficiency.

Since the focus of Ref. Lin  _et al._ , 2019 is to show the feasibility of the
quantum adiabatic doping for incommensurate filling, only one filling factor
at the particle-doped regime was studied. It is worth more systematic study
how the quantum adiabatic doping behaves from different fillings, and in
particular whether it remains efficient in the hole-doped regime should be
addressed.

In this work we carry out a systematic study on the quantum adiabatic doping
in a one-dimensional optical lattice for a broad range of filling factors,
with numerical simulations based on the time-dependent density matrix
renormalization group (DMRG) method. We find that the localization slowing
down is a generic problem for both commensurate and incommensurate fillings.
For the incommensurate case, the localization problem is more fundamental
because the localization persists in the thermodynamic limit. This problem
also causes slowing down for the commensurate case considering a finite-size
system with a localization length significantly smaller than the system size,
although it scales linearly with the system size in the thermodynamic limit.
For the particle doping, we show that the adiabatic preparation efficiency can
be enhanced by introducing atomic interaction for both commensurate and
incommensurate fillings. For the hole doping, we find that the localization is
much weaker, which we attribute to the large particle tunneling of the final
lattice. The quantum adiabatic hole doping is consequently more efficient than
the particle doping, and the efficiency can be further improved by including
strong atomic interaction. Besides starting from the band insulator, we also
consider the adiabatic particle doping starting from the Mott insulator at
half filling. Our numerical simulation shows that the quantum adiabatic doping
starting form the Mott insulator has better performance for certain fillings.
We expect these numerical results on the quantum adiabatic doping for a one-
dimensional optical lattice would also shed light on the two-dimensional
lattice.

##  II Theoretical Setup

The atomic quantum simulator of the Fermi-Hubbard model consists of two-
component Fermionic atoms confined in the periodical optical lattice. The
system is described by the Hamiltonian

$\displaystyle H$ $\displaystyle=$ $\displaystyle\int
d^{d}\mathbf{x}\psi^{\dagger}_{\sigma}\Big{(}-\frac{\hbar^{2}\nabla^{2}}{2M}+V(\mathbf{x})-\mu\Big{)}\psi_{\sigma}$
(1)
$\displaystyle+g\psi^{{\dagger}}_{\uparrow}\psi^{{\dagger}}_{\downarrow}\psi_{\downarrow}\psi_{\uparrow},$

where $\psi_{\sigma=\uparrow,\downarrow}({\bf x})$ is the quantum field
operator for corresponding pseudospin (hyperfine state) $\uparrow$ and
$\downarrow$ components, $M$ the atomic mass, $\mu$ the chemical potential,
$g$ the interaction strength between the two components, and $V(\mathbf{x})$
the optical lattice potential. We neglect the harmonic potential in our
calculation for simplicity. In this work we first consider a one-
dimensional($d=1$) band insulator as the initial state of the adiabatic
evolution. In optical lattice experiments, band insulators with low entropy
have been achieved Chiu  _et al._ (2018). The initial lattice potential reads
as

$V_{I}(x)=V\mathrm{cos}(2\pi x/\lambda).$ (2)

with $\lambda$ the lattice constant and $V$ the strength. Then we
adiabatically convert the lattice to another one with a different period

$V_{F}(x)=V^{\prime}\mathrm{cos}(2\pi x/\lambda^{\prime}).$ (3)

During the time evolution, the potential has a time dependent form

$V(x,t)=(1-s(t/T))V_{I}(x)+s(t/T)V_{F}(x),$ (4)

which is standard in the context of adiabatic algorithm with $s(t/T)$ giving
the path of the evolution and T the total evolution time Albash and Lidar
(2018). In the calculation, we focus at the parameter choice of
$V=V^{\prime}$, and examine the consequence of varying the overall lattice
strength. Assuming that the atomic loss during the evolution is negligible,
which holds when the evolution time $T$ is smaller than the lifetime of the
cold atom system, the filling factor of the final state is
$f=\lambda^{\prime}/\lambda$. By controlling the the ratio between the lattice
constants, a generic filling factor is accessible with the quantum adiabatic
doping.

With the incommensurate potential in Eq. (4), the quantum dynamics during the
adiabatic evolution cannot be described by a valid tight-binding model. We
thus have to take into account the continuous degrees of freedom of the
lattice. In the numerical simulation, we discretize the space as $x\to j\times
a$ , with $a$ the length of a grid, and $j$ the discrete index. This leads to
a Hamiltonian

$H=\sum_{j}\left\\{-t\left[c_{j\sigma}^{\dagger}c_{j+1,\sigma}+H.c.\right]+V_{j}c_{j\sigma}^{\dagger}c_{j\sigma}+Un_{j\uparrow}n_{j\downarrow}\right\\}$
(5)

with parameters $t=\hbar^{2}/(2Ma^{2})$, $V_{j}=V(ja)$ and $U=g/a$. In the
following calculation, we divide each spatial period of the initial lattice
into $20$ grids if not specified.

The linear path $s(t/T)=t/T$ is adopted for the adiabatic evolution. To
evaluate the performance of the adiabatic doping, we calculate the wave
function overlap

$\text{Overlap}=|\langle\Psi_{g}|\Psi_{T}\rangle|$ (6)

and the excitation energy

$\Delta
E=\langle\Psi(T)|H_{F}|\Psi(T)\rangle-\langle\Psi_{g}|H_{F}|\Psi_{g}\rangle,$
(7)

where $|\Psi_{g}\rangle$ is the ground state of the final Hamiltonian $H_{F}$,
and $|\Psi(T)\rangle$ is the final state of the quantum adiabatic evolution.
To confirm our finding is generic for finite-size system, we consider more
than one system sizes in following numerical study.

##  III Particle doping

###  III.1 Free fermion

Figure 1:  Quantum adiabatic particle doping of free fermions in one-
dimensional optical lattice with rational filling factors. The lines show the
dependence of the final state wave function overlap on total evolution time
$T$ with varying strengths of lattice potential. The filling factor after the
lattice conversion takes the rational values $f=3/4$ in (a),(b) and $f=2/3$ in
(c),(d). For each filling factor, we consider two different system sizes,
$L=36$ and $60$, with the results shown in (a),(c) and (b),(d), respectively.
The wave function overlap systematically increases with the adiabatic time
$T$. For weak lattice potential, the overlap approaches to $1$ with
sufficiently long evolution time and the doped final state can be efficiently
prepared. For stronger potential, the overlap remains nearly zero and the
adiabatic doping proposal fails for free fermions.

We first study the adiabatic particle doping of free fermions. This
corresponds to setting the filling factor $f$ larger than $1/2$. We simulate
the adiabatic evolution with rational filling factors $f=2/3,3/4$ and
calculate the wave function overlap of the final state with the ground state
(Eq. (6)). The dependence of the wave function overlap on the total evolution
time $T$ is shown in Fig. 1. The initial state is a one-dimensional band
insulator without interaction in optical lattice with $L$ periods. We consider
two different choices of system size, $L=36$, and $60$. For current and
following calculations, we use periodical boundary condition for non-
interacting fermions, while the open boundary condition is adopted for the
DMRG simulation of interacting fermions for numerical implementation
convenience. For weak lattice potential, say $V=1,2E_{R}$($E_{R}$ is the
single photon recoil energy), the overlap increases with $T$ and quickly
approaches to $1$ as we increase the total evolution time. This implies the
adiabatic preparation of the final state is efficient for the weak lattice
confinement. However, with strong lattice potential ($V=8E_{R}$, for example),
the final state overlap essentially remains at zero for all evolution time we
have simulated, which means the adiabatic doping is inefficient.

For the adiabatic doping with incommensurate lattice, the localization in the
intermediate dynamics leads to slowing down of the adiabatic state evolution.
This has been shown by the increase of inverse participation ratio(IPR) with
the lattice potential strength Lin  _et al._ (2019). In the commensurate case,
similarly, the breakdown of adiabatic state preparation under strong lattice
confinement indicates a localization-like problem also occurs during
evolution. To investigate the slowing down problem in the commensurate
lattice, we calculate the normalized participation ratio(NPR) Li  _et al._
(2017). The NPR of a single-particle eigenstate $\phi_{m}(x)$ is defined as

$\text{NPR}^{(m)}=\left[L\sum_{j}|\phi_{m}(j)|^{4}\right]^{-1},$ (8)

where $L$ is the system size, and $j$ labels the discrete spatial
coordinate(see Eq. (5)). This quantity remains finite for spatial extended
states but vanishes for localized states. For a one-dimensional localized
system with length $L$, it goes like $L^{-1}$. We calculate the single-
particle NPR of $L$ lowest eigenstates and average them to get

$\langle\text{NPR}\rangle=\frac{1}{L}\sum_{m=1}^{L}\text{NPR}^{(m)}.$ (9)

The averaged NPR is multiplied by $L$ to compensate the $L^{-1}$ scaling:

$M(L)=L\langle\text{NPR}\rangle.$ (10)

Therefore, for localized one-dimensional system, $M(L)$ is expected to be
nearly independent of $L$. In Fig. 2 we show the dependence of the quantity
$\text{log}_{10}M(L)$ on $L$ following the evolution path. The filling factor
is set to be $f=3/4$ and the overall potential strength is
$V=V^{\prime}=8\leavevmode\nobreak\ E_{R}$. It is evident that $M(L)$
increases quickly with $L$ at each point of the path. While this is consistent
with the well-known fact that commensurate lattice models do not have
localization in the thermodynamic limit, having a $\langle{\rm NPR}\rangle$
significantly smaller than 1 in a finite size system implies the coupling
between different modes is severely suppressed owing to the locality of the
Hamiltonian, which then causes the slowing down of the quantum adiabatic
evolution.

Figure 2:  The dependence of the quantity $\text{log}_{10}M(L)$ (Eq. (10)) on
the system size $L$ for the adiabatic particle doping of one-dimensional free
fermions. Here $s\in[0,1]$ is the path parameter of adiabatic evolution (Eq.
(4)). The system size $L$ takes values from $15$ to $597$ with a rational
filling factor $f=2/3$. In the calculation, the potential strength is
$V=V^{\prime}=8E_{R}$, for which the intermediate regime of evolution is
strongly localized. It is evident that this quantity increases quickly with
$L$ at each time instant $s$. For a one-dimensional localized system, $M(L)$
is expected to be nearly independent of $L$ due to the $L^{-1}$ scaling of
$\langle{\rm NPR}\rangle$. The increase of $\text{log}_{10}M(L)$ implies that
the localization becomes unstable as we increase the system size and tends to
disappear in the thermodynamic limit. Nonetheless, the significant deviation
of M(L) from L already implies the coupling between wave functions is rather
weak, which affects the quantum adiabatic doping.

###  III.2 Interacting fermion

The localization makes the minimal energy gap between the ground state and the
first excited state exponentially small, and therefore causes slowing down of
the adiabatic state preparation. In the study of many-body localization, it
has been established that interaction effect tends to destabilize localization
Basko  _et al._ (2006); Oganesyan and Huse (2007); Pal and Huse (2010); Kondov
_et al._ (2015); Schreiber  _et al._ (2015); Smith  _et al._ (2016); Choi  _et
al._ (2016), which would then improve the efficiency of the quantum adiabatic
doping. This has been shown to be efficient for an incommensurate lattice with
strong interaction in Ref. Lin  _et al._ , 2019. In this work, we consider
introducing a time-dependent interaction

$H_{\text{U}}=g(t)\int
dx\psi^{{\dagger}}_{\uparrow}\psi^{{\dagger}}_{\downarrow}\psi_{\downarrow}\psi_{\uparrow}$
(11)

with the time sequence shown in the inset of Fig. 3. The system is initially
prepared in a one-dimensional noninteracting band insulator with $14$ periods.
The interaction is turned on slowly and ramped to a constant. Then the initial
optical lattice is adiabatically converted to another spatial period, after
which the interaction is turned off slowly. The size of the final lattice is
$L^{\prime}=21$, which corresponds to a rational filling factor $f=2/3$. We
simulate the lattice conversion process using density matrix renormalization
group (DMRG) method. The total evolution time is $T=200\leavevmode\nobreak\
\hbar/E_{R}$. We use the second-order Trotterization White and Feiguin (2004)
of the evolution operator with $26000$ evolution steps. The wave function
overlap and the single-particle excitation energy at $t=T$ are shown as
functions of the interaction strength $g$ in Fig. 3. Here the overall
potential strength is $V=V^{\prime}=8E_{R}$, for which the non-interacting
quantum adiabatic doping is inefficient. As we increase the interaction
strength $g$ from $0$ to $0.2$ (in unit of $\lambda E_{R}$ with $\lambda$ the
spatial period of the initial lattice), the wave function overlap is improved
from $0.128$ to $0.484$, and the single-particle excitation energy $\Delta
E/N$ ($N$ is the number of particles) becomes much suppressed. We thus find
that in the commensurate case, the slowing down problem can still be solved by
introducing atomic interaction. An dramatically enhancement of the preparation
efficiency is achieved with a proper interaction strength. However, when we
further increase the interaction strength, we find that the overlap decreases
and $\Delta E/N$ increases slightly. In the dynamical simulation, the largest
bond dimension of the matrix product state(MPS) is $\chi=120$, with which
numerical convergence is reached (see Ref. Lin  _et al._ , 2019)

Figure 3:  Quantum adiabatic particle doping of one-dimensional interacting
fermions with rational filling. We simulate the adiabatic evolution with DMRG,
taking the Hamiltonian in Eq. (5). The inset shows the time sequence of the
evolution. The interaction is adiabatically tuned on for a one-dimensional
noninteracting band insulator, held constant to implement the lattice
conversion, and then adiabatically tuned off. We perform the calculation of
the lattice conversion process in the regime $t\in[0,T]$. In the DMRG
calculation, we take second-order Trotterization of the evolution operator,
with the evolution time $T=200\leavevmode\nobreak\ \hbar/E_{R}$ divided into
$26000$ Trotter steps. The filling factor is chosen as $f=2/3$ with $L=14$ and
$L^{\prime}=21$. The strength of the potential is
$V=V^{\prime}=8\leavevmode\nobreak\ E_{R}$. As we increase the interaction
strength from $g=0$ to $g=0.2\lambda E_{R}$, the final state wave function
overlap increases and the single-particle excitation energy $\Delta E/N$ ($N$
is particle number) becomes smaller. The localization problem is reduced and
the performance of the adiabatic doping is improved by introducing interaction
of proper strength.

We also consider performing quantum adiabatic doping starting from a Mott
insulating state. In experiments, such an initial state is accessible since
the low-temperature antiferromagnetic order has been observed at half filling
in two-dimensional optical lattices Mazurenko  _et al._ (2017). To compare the
performances of starting from the two different initial states of Mott and
band insulators, we simulate the adiabatic lattice conversion starting from
interacting band insulator with $L$ periods and Mott insulator with $2L$
periods. We choose five different fillings in the particle-doped regime, both
rational and irrational. Without loss of generality, the irrational filling
factor is set to be golden ratio, which is approximated by the Fibonacci
sequence as $8/13$ and $13/21$ in our finite-size calculation. Here we
consider several different system sizes, $L=8,9,12,13$. For the Mott
insulator, each period of the initial lattice is divided into $10$ grids,
while that for band insulator is divided into $20$ grids. This parameter
choice is chosen such that the final states of two procedures have the same
spatial periods and discrete grids. The total evolution time is
$T=200\leavevmode\nobreak\ \hbar/E_{R}$ for the smaller system of $L=8,9$ and
$T=200,300,400\leavevmode\nobreak\ \hbar/E_{R}$ with the corresponding numbers
of Trotter steps $26000$, $38000$, and $50000$ for the larger one of
$L=12,13$. During the evolution, the interaction strength is held constant of
$g=0.4\lambda E_{R}$ and the overall potential strength takes the value
$V=V^{\prime}=8\leavevmode\nobreak\ E_{R}$. The final state wave function
overlap is shown in Fig. 4 with the solid and dashed lines representing the
results for band insulator and Mott insulator, respectively. The overlap is
larger for band insulator at $f=3/5$ and $(\sqrt{5}-1)/2$ .

For $f=2/3$ and $3/4$, the adiabatic doping starting from Mott insulator has
better performance. The advantage is especially evident for $f=3/4$. As shown
in Fig. 4(c), the overlap reaches $0.93$ at $T=300\leavevmode\nobreak\
\hbar/E_{R}$ for Mott insulator while that for band insulator is $0.17$. At
$f=4/5$, the overlaps of the two protocols are both smaller than $1/2$ for all
adiabatic time we have simulated. As shown in Fig. 4(d), the overlap staring
from Mott insulator exceeds that of band insulator and reaches $0.31$ at
$T=400\leavevmode\nobreak\ \hbar/E_{R}$. We expect the overall performance can
be further improved by increasing the total evolution time.

Figure 4:  Performance of quantum adiabatic doping of one-dimensional
interacting fermions starting from band insulator and Mott insulator. We
simulate the adiabatic evolution using DMRG with several different system
sizes, $L=8,9,12,13$. For the Mott insulator, the corresponding number of
lattice spatial periods is $2L$ with the lattice constant $\lambda/2$. During
the lattice conversion, the interaction strength is held constant of
$g=0.4\lambda E_{R}$. The final state wave function overlap of these two
different protocols are compared in (a) for $L=8,9$ with the evolution time
$T=200\leavevmode\nobreak\ \hbar/E_{R}$ and (b) for $L=12,13$ with
$T=200,300,400\leavevmode\nobreak\ \hbar/E_{R}$. We choose $5$ different final
state fillings, both rational and irrational. In (c) and (d), we show the
dependence of wave function overlap on the adiabatic time $T$ with $f=3/4$ and
$4/5$, respectively. For $f=2/3$ and $3/4$, the adiabatic doping tarting from
the Mott insulator has better performance. The comparison is most dramatic for
$f=3/4$, as shown in (c). For $f=4/5$, the wave function overlaps of both
procedures are smaller than $1/2$ in the time regime we consider. As we
increase $T$, the protocol starting from a Mott insulator is more efficient
than that from a band insulator.

##  IV Hole doping

###  IV.1 Free fermion

In this section we study the adiabatic hole doping, which corresponds to the
filling factor smaller than $1/2$. We first simulate the adiabatic evolution
of free fermions starting from one-dimensional non-interacting band insulator
in an optical lattice with $L$ periods. Both the rational and irrational
fillings are considered. The rational filling factors are chosen as $f=1/3$
and $1/4$, while for the irrational case, we take
$f=\left[3-\sqrt{5}\right]/2$, which is approximated by Fibonacci sequence as
$f=F_{n-2}/F_{n}$. We choose four different values for the system size,
$L=36,60$ for the rational fillings and $L=34,55$ for the irrational case. The
dependence of the final state wave function overlap on the total evolution
time is shown in Fig. 5. The overlap systemically increases with evolution
time for all potential strengths we consider. It should be noticed that a wave
function overlap larger than $1/2$ can be reached for the larger system of
$L=55,60$ with $V=8E_{R}$ and smaller system of $L=34,36$ with $V=16E_{R}$ as
we increase the evolution time to $T=800\hbar/E_{R}$. The adiabatic doping
remains efficient in the hole-doped regime even for moderate lattice
potentials, for example with $V=8E_{R},16E_{R}$, in contrast to the particle-
doping, which is severely subjected to localization slowing down problem. For
sufficiently strong lattice potential, the slowing down problem still occurs
in the quantum adiabatic evolution of the hole-doped case.

Figure 5:  Quantum adiabatic hole doping of free fermions starting from one-
dimensional band insulator. The solid and dashed lines show the dependence of
final state wave function overlap on total evolution time $T$ for hole doping
and particle doping, respectively. For hole doping, the filling factor takes
rational values $1/3$ and$1/4$ in (a), (b) and (c), (d), respectively. The
system sizes are set to be $L=36$ in (a), (c) and $L=60$ in (b), (d). It also
takes irrational value $\left[3-\sqrt{5}\right]/2$, which is approximate by
the Fibonacci sequence as $f=34/89$ with $L=34$ in (e) and $55/144$ with
$L=55$ in (f). The non-interacting hole doping is efficient for weak lattice
confinement, say $V=4E_{R}$. That is, the wave function overlap quickly
approaches to $1$ for rational fillings we considered, and also the smaller
system of $L=34$ with irrational filling. Under strong lattice confinement
with $V=32E_{R}$, the overlap remains small and the quantum adiabatic doping
fails to prepare the final state. For the same strength, the adiabatic doping
is more efficient in the hole-doped regime. For the fillings in hole-doped
regime we consider, a overlap larger than $1/2$ can be reached at
$T=800\hbar/E_{R}$ for $L=34,36$ with $V=16E_{R}$ and $L=55,60$ with
$V=8E_{R}$, while particle doping is inefficient with $V=8E_{R}$ (shown by
dashed lines with the corresponding filling factors $1-f$).

To investigate the localization in the hole-doped regime, we calculate the
inverse participation ratio(IPR), which tends to vanish in the extended system
and remains finite in the localized system. The IPR is averaged over the $L$
lowest lying single-particle eigenstates following the evolution path. The
results are shown in Fig. 6. The irrational filling is set to be
$f=\left[3-\sqrt{5}\right]/2$, which is approximated by $55/144$ with $L=55$.
The averaged IPR systematically increases with the potential strength, which
means the breakdown of adiabatic preparation corresponds to the atom
localization. The results are compared to that in the particle-doped regime
with $f=\left[\sqrt{5}-1\right]/2$ (approximated by $55/89$) and $V=8E_{R}$.
For the most part of Hamiltonian evolution path, the IPR of hole doping is
smaller than that of particle doping, especially for $s>0.5$ where the
strength of final potential becomes dominant. This is because in comparison
with the particle-doping case, the final lattice constant is smaller for hole
doping, for which quantum tunneling is stronger, suppressing the atom
localization. The quantum adiabatic doping protocol is evidently more
efficient in the case of hole-doping than the particle doping.

Figure 6:  Averaged inverse participation ratio(IPR) along the evolution path
with $V=8,16,32E_{R}$ in the hole-doped regime and $V=8E_{R}$ in the particle-
doped regime. The filling factor is set to be $f=\left[3-\sqrt{5}\right]/2$
and $\left[\sqrt{5}-1\right]/2$ and approximated by Fibonacci sequence as
$F_{n-2}/F_{n}$ and $F_{n-2}/F_{n-1}$ with the initial lattice size $L=55$,
respectively. The averaged IPR increases with potential strength, which
implies the localization becomes stronger. Compared to particle doping, the
IPR for hole doping is smaller for most part of the evolution path, which
explains the high preparation efficiency in this regime.

###  IV.2 Interacting fermion

Figure 7:  Performance of adiabatic hole doping of interacting fermions in
one-dimensional lattice for both rational ($f=2/5$) and irrational
fillings($\left[3-\sqrt{5}\right]/2$). We simulate the dynamics using DMRG
with the same parameter choice as in Fig 3. The wave function overlap and
single particle excitation energy are shown with solid lines in (a),(b) for
rational filling with $L=12$ and (c),(d) for irrational filling, which is
approximated by $L/L^{\prime}=13/34$. Here we consider the total adiabatic
time $T=50\hbar/E_{R}$ and $200\hbar/E_{R}$. The overall potential strength is
$V=8E_{R}$ for both particle and hole doping. It is found that the preparation
efficiency is improved by introducing strong interaction. As we increase the
interaction strength, the wave function increases and approaches to $1$ for
$T=200\hbar/E_{R}$, and the excitation energy is reduced accordingly. Compared
to particle doping (shown by dashed lines with filling factor $f=13/21$), the
hole doping shows a better preparation efficiency for this interacting case,
as is true for free fermions as well.

We further consider adiabatic hole doping of interacting fermions, and
simulate the evolution process starting from a one-dimensional interacting
band insulator using DMRG. We choose $f=2/5$ with $L=12$ for the rational
filling and $f=\left[3-\sqrt{5}\right]/2$ for the irrational one, which is
approximated by $L/L^{\prime}=13/34$ in the calculation. The dependence of the
final state wave function overlap and single particle excitation energy on the
interaction strength is shown in Fig. 7 (solid lines) with an overall
potential strength $V=8E_{R}$. The DMRG calculation shares the same time
sequence and parameter choice with that in the particle-doped regime. Two
choices of total evolution time, $T=50\hbar/E_{R}$ and $200\hbar/E_{R}$, are
considered with $7000$ and $26000$ evolution steps, respectively. The
preparation efficiency is improved by introducing strong atomic interaction.
In our numerical results, we find systematic increase of wave function overlap
(Eq. (6)) and reduction of excitation energy with increasing interaction
strength. By increasing interaction strength $g$ from $0$ to $0.6$ (in unit of
$\lambda E_{R}$), the wave function overlap ultimately reaches $0.96$ for
$f=13/34$ and $0.965$ for $f=2/5$. We compare the results with that of
particle doping (dashed line) for the same interaction and confinement
potential strengths. As shown in Fig. 7 (c), the wave function overlap for
$T=50\hbar/E_{R}$ in the hole-doped regime is even large than that of particle
doping for $T=200\hbar/E_{R}$. The hole doping is evidently more efficient
than the particle doping in the interacting regime, as is true for non-
interacting case as well.

##  V Conclusion

To conclude, the doped quantum phase of atomic Fermi-Hubbard model with low
thermal entropy can be prepared by adiabatically converting two optical
lattices with different spatial periods. In this adiabatic doping proposal, an
arbitrary filling fraction can be achieved by choosing the lattice constant
ratio of the initial and final lattices. In this work, we consider the quantum
adiabatic doping starting from band insulator in one-dimensional lattice and
systemically study the proposal for a broad range of filing fractions from
particle-doping to hole-doping, including both rational and irrational cases.
It is found that the atom localization, which is a fundamental problem in
incommensurate lattice, also prevents efficient adiabatic doping for
commensurate filling at strong lattice confinement. Through DMRG simulation,
we show that the localization is suppressed by introducing strong atomic
interaction and the state preparation efficiency is consequently improved.
Compared to particle doping, the adiabatic hole doping is more efficient in
both free and interacting regime. We also consider the adiabatic doping
starting from the Mott insulator state at half filling. From the results of
DMRG calculation, we find that this protocol has significantly higher
preparation efficiency than the adiabatic doping starting from band insulator
for a certain range of filling factors.

Although the numerical simulation of the quantum adiabatic doping process is
restricted to one dimension in this work due to numerical cost, we anticipate
a better performance in two dimensions for the localization physics is largely
expected to be weaker in higher dimensions. Our systematic study of different
fillings and interaction effects should also shed light on the quantum
adiabatic doping in two dimensions.

Acknowledgement. This work is supported by National Key R&D Program of China
(under grant no. 2018YFA0306502), National Natural Science Foundation of China
(Grants No. 11934002 and 11774067), National Program on Key Basic Research
Project of China (Grant No. 2017YFA0304204), and Shanghai Municipal Science
and Technology Major Project (Grant No. 2019SHZDZX01). Calculations were
performed based on the ITensor Library ITe .

## References

  * Lin  _et al._ (2019) J. Lin, J. Nan, Y. Luo, X.-C. Yao, and X. Li, Physical Review Letters 123, 233603 (2019).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Reviews of modern physics 80, 885 (2008).
  * Li and Liu (2016) X. Li and W. V. Liu, Reports on Progress in Physics 79, 116401 (2016).
  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Advances in Physics 56, 243 (2007).
  * Bloch (2018) I. Bloch, Nature Physics 14, 1159 (2018).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Hofstetter  _et al._ (2002) W. Hofstetter, J. I. Cirac, P. Zoller, E. Demler, and M. Lukin, Physical review letters 89, 220407 (2002).
  * Dutta  _et al._ (2015) O. Dutta, M. Gajda, P. Hauke, M. Lewenstein, D.-S. Lühmann, B. A. Malomed, T. Sowiński, and J. Zakrzewski, Reports on Progress in Physics 78, 066001 (2015).
  * Lee  _et al._ (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Reviews of modern physics 78, 17 (2006).
  * Metzner  _et al._ (2012) W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Reviews of Modern Physics 84, 299 (2012).
  * Jiang and Devereaux (2019) H.-C. Jiang and T. P. Devereaux, Science 365, 1424 (2019).
  * Qin  _et al._ (2020) M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 10, 031016 (2020).
  * Greif  _et al._ (2013) D. Greif, T. Uehlinger, G. Jotzu, L. Tarruell, and T. Esslinger, Science 340, 1307 (2013).
  * Hart  _et al._ (2015) R. A. Hart, P. M. Duarte, T.-L. Yang, X. Liu, T. Paiva, E. Khatami, R. T. Scalettar, N. Trivedi, D. A. Huse, and R. G. Hulet, Nature 519, 211 (2015).
  * Mazurenko  _et al._ (2017) A. Mazurenko, C. S. Chiu, G. Ji, M. F. Parsons, M. Kanász-Nagy, R. Schmidt, F. Grusdt, E. Demler, D. Greif, and M. Greiner, Nature 545, 462 (2017).
  * Trebst  _et al._ (2006) S. Trebst, U. Schollwöck, M. Troyer, and P. Zoller, Phys. Rev. Lett. 96, 250402 (2006).
  * Sørensen  _et al._ (2010) A. S. Sørensen, E. Altman, M. Gullans, J. Porto, M. D. Lukin, and E. Demler, Physical Review A 81, 061603 (2010).
  * Lubasch  _et al._ (2011) M. Lubasch, V. Murg, U. Schneider, J. I. Cirac, and M.-C. Banuls, Physical review letters 107, 165301 (2011).
  * Zhang and Duan (2013) Z. Zhang and L.-M. Duan, Phys. Rev. Lett. 111, 180401 (2013).
  * Chiu  _et al._ (2018) C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, and M. Greiner, Phys. Rev. Lett. 120, 243201 (2018).
  * Sun  _et al._ (2020) H. Sun, B. Yang, H.-Y. Wang, Z.-Y. Zhou, G.-X. Su, H.-N. Dai, Z.-S. Yuan, and J.-W. Pan, arXiv preprint arXiv:2009.01426 (2020).
  * Basko  _et al._ (2006) D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annals of physics 321, 1126 (2006).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Physical review b 75, 155111 (2007).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Physical Review B 82 (2010).
  * Kondov  _et al._ (2015) S. S. Kondov, W. R. McGehee, W. Xu, and B. DeMarco, Phys. Rev. Lett. 114, 083002 (2015).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, Reviews of Modern Physics 90, 015002 (2018).
  * Li  _et al._ (2017) X. Li, X. Li, and S. D. Sarma, Physical Review B 96, 085119 (2017).
  * Schreiber  _et al._ (2015) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, Science 349, 842 (2015).
  * Smith  _et al._ (2016) J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nature Physics 12, 907 (2016).
  * Choi  _et al._ (2016) J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, Science 352, 1547 (2016).
  * White and Feiguin (2004) S. R. White and A. E. Feiguin, Physical review letters 93, 076401 (2004).
  * (32) ITensor Library (version 2.0.11) http://itensor.org .