Quantum Adiabatic Doping with Incommensurate Optical Lattices

TL;DR

This paper proposes a method for doping Fermi-Hubbard models in optical lattices via adiabatic evolution, addressing localization challenges and suggesting interaction-induced delocalization to improve efficiency in low-temperature quantum simulations.

Contribution

It introduces a theoretical protocol for adiabatic doping in incommensurate optical lattices, overcoming localization issues with many-body delocalization techniques.

Findings

Localization hinders adiabatic doping in strong lattice regimes.

Interaction-induced delocalization can mitigate slowing down in 1D.

The protocol is potentially effective in 2D where localization is less stable.

Abstract

Quantum simulations of Fermi-Hubbard models have been attracting considerable efforts in the optical lattice research, with the ultracold anti-ferromagnetic atomic phase reached at half filling in recent years. An unresolved issue is to dope the system while maintaining the low thermal entropy. Here we propose to achieve the low temperature phase of the doped Fermi-Hubbard model using incommensurate optical lattices through adiabatic quantum evolution. In this theoretical proposal, we find that one major problem about the adiabatic doping that shows up is atomic localization in the incommensurate lattice, potentially causing exponential slowing down of the adiabatic procedure. We study both one- and two-dimensional incommensurate optical lattices, and find that the localization prevents efficient adiabatic doping in the strong lattice regime for both cases. With density matrix renormalization group calculation, we further show that the slowing down problem in one dimension can be circumvented by considering interaction induced many-body delocalization, which is experimentally feasible using Feshbach resonance techniques. This protocol is expected to be efficient as well in two dimensions where the localization phenomenon is less stable.

Figures

Equations

$$H=\int d^d\mathbf{x}{\psi }_{\sigma }^{†}\left[-\frac{{\hslash }^2{\vec{\nabla }}^2}{2M}+V(\mathbf{x})-\mu \right]{\psi }_{\sigma }+g{\psi }_{\uparrow }^{†}{\psi }_{\downarrow }^{†}{\psi }_{\downarrow }{\psi }_{\uparrow },$$

$$V_I(\mathbf{x})=V{\sum }_{\alpha }\cos (2\pi x_{\alpha }/\lambda ),$$

$$V_F(\mathbf{x})=V^{\prime }{\sum }_{\alpha }\cos (2\pi x_{\alpha }/{\lambda }^{\prime }),$$

$$V(\mathbf{x},t)=[1-s(t/T)]V_I(\mathbf{x})+s(t/T)V_F(\mathbf{x}),$$

$$H={\sum }_j\left\{-t\left[c_{j\sigma }^{†}{c}_{j+1,\sigma }+H.c.\right]+V_j{c}_{j\sigma }^{†}{c}_{j\sigma }+U{n}_{j\uparrow }{n}_{j\downarrow }\right\},$$

$$\mathrm{Overlap}=|⟨{\Psi }_g|{\Psi }_{\mathrm{final}}⟩|,$$

$$H(t)=K+V(t),$$

$$|\Psi (t)⟩={\psi }_1^{†}(t){\psi }_2^{†}(t)\dots {\psi }_N^{†}(t)|0⟩,$$

$${\psi }_j^{†}=\sum _x{U}_{xj}(t)c_x^{†}.$$

$$i{\partial }_t{U}_{xj}(t)=\sum _{x^{\prime }}\left[{\mathcal{K}}_{x{x}^{\prime }}+{\mathcal{V}}_{x{x}^{\prime }}(t)\right]U_{x^{\prime }j}(t).$$

$$U(t+\delta t)=e^{-i\mathcal{K}\delta t/2}{e}^{-i\mathcal{V}(t+\delta t/2)\delta t}{e}^{-i\mathcal{K}\delta t/2}U(t)+\mathcal{O}(\delta t^3).$$

$$E(t)=\sum _{j=1}^N\sum _{x{x}^{\prime }}\left[{\mathcal{T}}_{x{x}^{\prime }}+{\mathcal{V}}_{x{x}^{\prime }}(t)\right]U_{xj}^{\ast }(t)U_{x^{\prime }j}(t).$$

$$\mathrm{Overlap}=|⟨{\Psi }_g|\Psi (t=T)⟩|.$$

$$\mathrm{Overlap}=|\mathrm{Det}[U^{†}{U}^{(g)}]|.$$

Full text

Introduction.— Quantum simulation with ultracold atoms confined in optical lattices have attracted tremendous efforts in the last decade Bloch (2018). One ultimate goal is to study quantum many-body physics at low temperature Jaksch et al. (1998); Hofstetter et al. (2002); Lewenstein et al. (2007); Bloch et al. (2008); Dutta et al. (2015); Li and Liu (2016) whose simulation on classical computers is subjected to unreachable computational complexity. In particular, quantum simulation of the low-temperature phase diagram of the Fermi-Hubbard model Hofstetter et al. (2002); Esslinger (2010) is of great experimental interest, as it would help understand the strongly correlated physics of direct relevance to making high-temperature superconductors 200 (2007). Owing to the high controllability of the microscopic degrees of freedom in the optical lattice emulator, the question of how the strongly correlated physics emerges can be addressed in an unarguable manner.

In theory, it is now well-known that Fermi-Hubbard like models could support various exotic quantum phases even in the weakly interacting regime. The combination of Fermi-Hubbard model with synthetic gauge fields leads to quantum Hall like topological states Osterloh et al. (2005); Ruseckas et al. (2005); Liu et al. (2009); Dalibard et al. (2011); its extension to multi-orbital setting gives rise to rich Fermi-surface nesting effects causing plentiful spontaneous symmetry breaking orders Hauke et al. (2011); Liu et al. (2014); Li and Liu (2016); its incorporation of long rang interactions supports unconventional density waves Li and Sarma (2015). These theoretical results obtained at weak interaction indicate the strongly interacting regime which is beyond classical simulation capability may contain even more fascinating physics Raghu and Kivelson (2011). To study such exotic physics demands the experimental optical lattices to enter the extremely low temperature, or more precisely the low entropy regime.

In the last few years, spectacular progress has been made in optical lattice Fermi-Hubbard emulator—the low-temperature antiferromagnetic phase at half filling has been reached Greif et al. (2013); Hart et al. (2015); Mazurenko et al. (2017) with quantum microscope techniques Haller et al. (2015); Cheuk et al. (2015); Parsons et al. (2015); Edge et al. (2015); Omran et al. (2015); Greif et al. (2016); Cheuk et al. (2016a); Parsons et al. (2016); Boll et al. (2016); Cheuk et al. (2016b); Brown et al. . With the experimental developments, Fermi surface nesting related many-body effects now become accessible. However, reaching the low entropy region of the doped Fermi-Hubbard model, in order to emulate strongly correlated electronic quantum physics such as the d-wave superconductivity, is experimentally challenging .

In this letter, we propose to use incommensurate optical lattices for quantum adiabatic doping. Preparing an initial band insulator state in a periodic optical lattice, the quantum state is adiabatically converted into another lattice having a different period. For a rational filling factor, adiabatic quantum simulation of the Fermi-Hubbard model can be achieved by considering a commensurate superlattice Trebst et al. (2006); Lubasch et al. (2011). To achieve a generic filling in the final optical lattice, the superposed lattice in the intermediate adiabatic process is necessarily incommensurate. We study both one- and two-dimensional incommensurate optical lattices, and find that a localization problem occurs when the lattice potential is strong. This problem leads to an exponentially small energy gap, and consequently prevents efficient adiabatic doping for both cases. With density matrix renormalization group calculation, we further show that the slowing down problem in one dimension can be circumvented by considering interaction induced many-body delocalization, which is experimentally feasible using Feshbach resonance techniques. This protocol is expected to be efficient as well in two dimensions where the localization phenomenon is less stable.

Quantum adiabatic doping.— A Fermi-Hubbard optical lattice emulator contains a two-component confined Fermi gas described by the Hamiltonian

[TABLE]

with $M$ the atomic mass, $\mu$ the chemical potential, $V({\bf x})$ the confining optical lattice, and $g$ the interaction strength between the two components. In this work, we focus on one- and two-dimensional cases, with $d=1$, and $2$, respectively. For simplicity, we have neglected harmonic trap potential in the calculation. In experiments, the trap potential should be compensated to keep evaporative cooling in the lattice Mathy et al. (2012); Duarte et al. (2015). The initial state we consider is a band insulator in a periodic lattice

[TABLE]

with $x_{\alpha}=x$ for a one-dimensional lattice, and $\alpha$ labelling the two orthogonal directions for the two dimensional case. A relevant energy unit is the single-photon recoil energy $E_{R}=\frac{\pi^{2}\hbar^{2}}{2M\lambda^{2}}$. With an optical lattice, a band insulator with low thermal entropy is experimentally accessible Chiu et al. (2018) for the energy gap of the band insulator ground state can be made as large as tens of kHz Bloch et al. (2008). We then adiabatically convert the system into another lattice having a different period,

[TABLE]

with a lattice constant $\lambda^{\prime}=\nu\lambda$ (see Fig. 1(a) for an illustration). During the adiabatic evolution, the lattice potential is time-dependent,

[TABLE]

with $s(t/T)=t/T$ a standard form of schedule as in the quantum adiabatic algorithm study Albash and Lidar (2018). Here $T$ is the total adiabatic time. The filling factor (the averaged particle number per degree of freedom) of the final lattice is $\nu^{d}$, assuming that the time of the adiabatic evolution is much shorter than the atom-loss timescale. Choosing different lattice constant ratio between the initial and final lattices, an arbitrary filling factor can be achieved. The particle and hole doping are achieved by setting $\nu^{d}$ greater and smaller than $1/2$, respectively.

Quantum adiabatic evolution of non-interacting fermions.— To demonstrate whether this quantum adiabatic doping works and estimate the required adiabatic time, we first investigate non-interacting fermions in one dimension (see Supplementary Material for the method). In solving the adiabatic evolution (see Eq. (4)), the tight-binding approximation is not applicable considering the intermediate time region. We thus have to take into account the continuous degrees of freedom of the lattice. In our calculation, the space coordinate is discretized as $x\to j\times a$, with $a$ the grid spacing. A periodic boundary condition is adopted to minimize finite-size effects. Without loss of generality, we set $\nu$ to be the golden ratio $[\sqrt{5}-1]/2$. In the finite size calculation with $L$ number of periods $\lambda_{1}$, $\nu$ is approximated using Fibonacci sequence $\{F_{n}\}$, as $\nu\approx F_{n-1}/F_{n}$, and $L=F_{n-1}$. The resultant lattice model is

[TABLE]

with $t=\frac{\hbar^{2}}{2Ma^{2}}$, $V_{j}=V(ja)$, and $U=g/a$. The interaction term $U$ is [math] for free fermions. In the numerical calculation, we divide each period of the original lattice potential evenly into $20$ grids. We expect the physics presented below to be generic against different choices of irrational filling or boundary conditions (Supplementary Materials).

The results are shown in Fig. 1. The wave function overlap between the final state of the adiabatic evolution ($|\Psi_{\rm final}\rangle$) and the ground state of the final Hamiltonian ($|\Psi_{g}\rangle$) is defined as

[TABLE]

with the result shown in Fig. 1(b). Fig. 1(c) shows the excitation energy in the final state as compared to the ground state. It is evident that the adiabatic quantum evolution is efficient in preparing the final ground state when the lattice is not too deep, say $V^{\prime}/E_{R}=1,2,4$. For a deep lattice, e.g. with $V^{\prime}/E_{R}=16$, the adiabatic evolution is found to be no longer efficient—the final state wave function overlap is approximately [math] and the excitation energy is significant. In Fig. 1(d), we show the inverse participation ratio, which is finite (vanishing) in the localized (extended) phase Iyer et al. (2013); Li et al. (2017). The ‘dashed’ lines in this plot correspond to the parameter-path of the adiabatic evolution. We find that the path enters into the localized regime for the inefficient adiabatic evolution at deeper lattice depths. This shows that the breakdown of the adiabatic preparation corresponds to the atom localization in the intermediate dynamics. In presence of localization, the level repulsion or the avoided crossing disappears due to the emergent local integrals of motion Serbyn et al. (2013); Huse et al. (2014); Chandran et al. (2015); Ros et al. (2015), and the minimal energy gap between the ground and first excited states becomes exponentially small. The localization thus causes the breakdown of the proposed adiabatic doping process. Also worth mentioning here is that the localization problem cannot be resolved by upgrading the linear ramp function $t/T$ in Eq. (4) to nonlinear ones.

We further consider the two-dimensional case, which is also of great experimental interest Greif et al. (2013); Hart et al. (2015); Mazurenko et al. (2017). In Fig. 2, we show the results of wave function overlap and excitation energy following the quantum adiabatic procedure using a two-dimensional incommensurate optical lattice (see Supplementary Material for the method). We still choose $\nu$ to be the golden ratio, so the eventual filling of the two dimensional lattice is $\nu^{2}=[\sqrt{5}-1]^{2}/4$. Like in the one-dimensional case, for a lattice not too deep, say with $V/E_{R}\leq 2$, the final state wave function overlap is close to $1$ and the excitation energy is tiny with a reasonable choice of adiabatic time, and the adiabatic doping is thus efficient. With a deeper lattice, the doping process becomes less efficient and eventually fails to prepare the doped ground state due to localization (see IPR in Supplementary Material). The localization problem is thus generic to the protocol of quantum adiabatic doping with incommensurate lattices, regardless of the dimensionality.

Atomic interaction comes to rescue the adiabatic doping.— To rescue the adiabatic doping against localization, one can perform Hamiltonian path optimization through the standard quantum control methods such as GRAPE Khaneja et al. (2005); Machnes et al. (2011); De Fouquieres et al. (2011), CRAB Caneva et al. (2011); Doria et al. (2011), or Krotov algorithms Eitan et al. (2011); Reich et al. (2012), which is expected to be helpful for experimental implementation of the quantum adiabatic doping. But the problem with that approach is that it is not generically applicable and may require optimization case-by-case for different lattice depths.

A generic approach to solve the localization slowing down problem is offered by the theoretical study of many-body localization—interaction effects tend to generically destabilize localization Basko et al. (2006); Schreiber et al. (2015); Kondov et al. (2015); Smith et al. (2016); yoon Choi1 et al. (2016). We thus introduce an auxiliary Hamiltonian beyond the standard adiabatic quantum computing algorithm. It has been shown in the context of quantum algorithms that this approach could lead to exponential speedup Van Dam et al. (2001); Aharonov et al. (2008); Yu et al. (2018). The auxiliary Hamiltonian we apply here to optical lattice system is the atomic interaction $\textstyle H_{\rm AU}=g(t)\int d^{d}{\bf x}\psi_{\uparrow}^{\dagger}\psi_{\downarrow}^{\dagger}\psi_{\downarrow}\psi_{\uparrow}.$ The time sequence we consider in this work is shown in the inset of Fig. 3. We first prepare a noninteracting band insulator in the initial lattice. Then we turn on the interaction via Feshbach resonance techniques, perform the adiabatic lattice conversion, and then switch off the interaction. The interaction strength during the adiabatic lattice conversion has to be large enough to disable the localization. For experiments targeting doped Fermi Hubbard model ground state with a particular interaction, one can eventually adiabatically tune the final interaction to that strength. We remark here that this work here is to show interaction effects could enable efficient adiabatic doping, rather than to optimize the time sequence.

Accelerating the one-dimensional adiabatic doping.— To explicitly demonstrate the feasibility of the adiabatic doping with interaction, we simulate the quantum dynamics of the one dimensional case with density matrix renormalization group (DMRG). Our DMRG program is developed according to the method in Ref. Schollwöck (2011) using matrix product stats. For numerical implementation convenience, we choose an open boundary condition in the DMRG calculation. Fig. 3 shows the the simulated results with $V/E_{R}=V^{\prime}/E_{R}=8$, corresponding to one parameter choice in Fig. 1 where localization occurs in absence of interaction. In the DMRG calculation, each period of the lattice potential is still divided into $20$ pieces as in the noninteracting case, but the number of periods simulated is reduced to $L=13$, due to expensive numerical cost. The dynamics is calculated using the second-order Trotterization with 26000 evolution steps for the adiabatic time of $T=200\hbar/E_{R}$. We choose the bond dimension to be $150$, for which the results of final state wave function overlap and excitation energy have converged (see Supplementary Material). In the regime where localization causes the breakdown of the adiabatic doping procedure, we find that interactions could help significantly enhance the final state wave function overlap and reduce the excitation energy. The performance of the adiabatic procedure can be further improved by increasing the time of adiabatic evolution. After the many-body localization in the adiabatic procedure is suppressed, the required adiabatic time to reach certain level of excitation energy is expected to scale polynomially with the system size for the many-body level repulsion of a thermal quantum system.

Discussion of the two-dimensional case.— For the two-dimensional case, although we could not afford to carry out the calculation due to the computational challenge—this challenge is the precise reason why quantum simulation of Fermi-Hubbard model is needed, we anticipate the interaction also solves the localization problem in two dimensions. In the study of many-body localization Basko et al. (2006), it is expected that the localization phenomenon is less stable in higher dimensions according to the thermal bubble argument yoon Choi1 et al. (2016); De Roeck and Huveneers (2017); Potirniche et al. (2018); Wahl et al. (2019). The protocol of adiabatic doping combined with auxiliary interaction is thus expected to hold as well for a two-dimensional lattice.

Experimental timescales.— Commonly used atoms to simulate Fermi-Hubbard model in optical lattice experiments are 40K and 6Li. Considering these atoms confined in an optical lattice formed by laser beams of wavelength $1064$ nm, the recoil energies are $E_{R}/\hbar\approx 2\pi\times 4$ kHz and $2\pi\times 29$ kHz, respectively. The required adiabatic time for the excitation energy per atom to drop down to one percent of single-particle tunneling is within $10$ms according to our calculation.

Conclusion.— To conclude, in quantum simulations of Fermi-Hubbard model with optical lattices, the adiabatic doping can be achieved by using incommensurate lattices. For noninteracting fermions, the localization problem makes the adiabatic doping inefficient in both one- and two-dimensional lattices in the strong potential region. As a generic recipe to solve this localization problem, we show that this problem can be solved in a generic manner by introducing strong atomic interactions through Feshbach resonance techniques, for interaction mediated atomic scattering causes many-body delocalization. This is confirmed by DMRG calculation for the one-dimensional system and is expected to hold for the two-dimensional case as well. These theoretical results may improve optical lattice quantum simulations of Fermi-Hubbard models in the doped regime of great importance to modeling high temperature superconductors or strongly correlated electrons in general.

Acknowledgement.— We acknowledge helpful discussion with Peter Zoller, Immanuel Bloch, Ana Maria Rey, Andrew Daley, Youjin Deng, and Markus Greiner. This work is supported by National Program on Key Basic Research Project of China under Grant No. 2017YFA0304204, 2018YFA0306501, National Natural Science Foundation of China under Grants No. 11774067, 11934002, 11874340. Calculations were performed based on the ITensor Library ITe . X.L would like to thank Department of Physics at Harvard University for hospitality during the completion of this work.

*†*Jian Lin and Jue Nan contributed equally to this work.

S-1 Numerical method to simulate a many-body state of free fermions

In this section, the details of our method in simulating non-interacting fermions are provided. The method is described by assuming a lattice Hamiltonian, with no loss of generality, because the Hamiltonian in continuum can be mapped to a lattice Hamiltonian through discretization—the problem of our interest here has no fermion doubling issue.

Considering a lattice model, with its degrees of freedom described by creation/annihilation operators ($c_{x}$/$c_{x}^{\dagger}$)—the subscript $x$ is the lattice-site index. For non-interacting fermions, the Hamiltonian takes a form of

[TABLE]

where $K=\sum_{xx^{\prime}}{\cal K}_{xx^{\prime}}c_{x}^{\dagger}c_{x^{\prime}}$ represents the time independent kinetic tunnelling, and $V(t)=\sum_{x}{\cal V}_{xx^{\prime}}(t)c_{x}^{\dagger}c_{x^{\prime}}$ the time dependent potential. The matrix ${\cal V}$ is diagonal. The time-evolving many-body state is described as

[TABLE]

where $N$ is the total particle number of fermions, and the operators $\psi_{1,2,\ldots N}$ are related to lattice operators by

[TABLE]

The $U$ matrix contains the single-particle wave function of each occupied orbital. From the many-body schrödinger equation, $i\partial_{t}|\Psi(t)\rangle=H(t)|\Psi(t)\rangle$, the $U$ matrix satisfies the single-particle schrödinger equation,

[TABLE]

In simulating the dynamics, we solve the time-dependent problem numerically by second-order Trotterization,

[TABLE]

Note that the matrix ${\cal V}$ is diagonal, so the multiplication of $e^{-i{\cal V}\delta t}$ on a vector can be implemented as multiplication of each element by a phase, which is numerically very efficient. The matrix ${\cal K}$ is non-diagonal, but can be diagonalized by a Fourier transformation. Then the multiplication of $e^{-i{\cal K}\delta t}$ can be efficiently carried out by using Fast-Fourier-Transformation. Regarding the initial condition, at time $t=0$, the initial state in our adiabatic doping procedure is the instantaneous ground state of the Hamiltonian at time $t=0$.

The energy of the dynamical state is defined with respect to the instantaneous Hamiltonian, which can be calculated according to

[TABLE]

The excitation energy of the adiabatic procedure is defined to be energy of the final state subtracted by the ground state energy of the final Hamiltonian. We also calculate the wave function overlap between the final state of the dynamical evolution and the ground state of the final Hamiltonian, defined as

[TABLE]

With the wave function of the occupied orbitals for $|\Psi_{g}\rangle$ stored as $U_{xj}^{(g)}$, the wave function overlap reads as

[TABLE]

For the spinful case, the two spin components can be further incorporated in a straightforward way.

S-2 Evolution of density profile in the quantum adiabatic doping

In our simulation of free fermions in the quantum adiabatic doping, besides the final state fidelity and excitation energy, we also calculate the evolution of the density profile. The results are shown in Fig. S1. When the lattice is shallow ($V=E_{R}$), the density profile follows the instantaneous ground state evolution, and the confined atoms slowly migrate from original lattice sites to the final lattice sites, and eventually becomes evenly distributed in the final lattice. In contrast, for a strong lattice ($V=8E_{R}$), the quantum tunneling from the original sites to the final sites is not efficient enough to make an even density profile when the adiabatic doping finishes.

S-3 Robustness against different doping fraction and boundary conditions

In optical lattice experiments, we typically do not have periodic boundary condition as used for the results of free fermions in the main text. To confirm our proposing quantum adiabatic doping is adaptable, we check its robustness against the choice of boundary conditions. From the results shown in Fig. S2, it is evident that for free fermions the quantum adiabatic doping remains efficient when the lattice is shallow, and becomes inefficient in a strong lattice. The behavior of the fidelity is qualitatively the same as the case of periodic boundary condition as shown in Fig. 1 (see main text). We thus expect the physics presented in this work to be largely independent of the choice of boundary conditions.

To confirm the quantum adiabatic doping works for generic fillings, we provide the results for the filling factor $\nu=\sqrt{2}$, which can be asymptotically approached by using the Pell’s number sequence. Despite the quantitative difference form the the choice of golden ratio, the generic feature—the efficiency (breakdown) at weak (strong) lattice—remains (see Fig. S3). We expect the physics in our quantum adiabatic doping to be generic for different filing factors.

S-4 Inverse participation ratio for the 2d lattice

Here we provide results of inverse participation ratio of the 2d lattice in the quantum adiabatic doping procedure (see Fig. S4). It is evident that the breakdown of the quantum adiabatic doping for the 2d case at strong lattice is also due to localization physics as in the 1d case. We thus expect the interaction restores the efficiency of the quantum adiabatic doping in two dimensions. In the study of many-body localization Basko et al. (2006), it is expected that the localization phenomenon is less stable in higher dimensions according to the thermal bubble argument yoon Choi1 et al. (2016); De Roeck and Huveneers (2017); Potirniche et al. (2018); Wahl et al. (2019).

S-5 Details of density-matrix-renormalization-group simulation

In this section, we provide the details of simulating interacting fermions with the density-matrix-renormalization-group (DMRG) method. In the adiabatic evolution, we represent the dynamical state with a matrix product state of a finite bond dimension $D$, and use the standard two-site procedure in simulating the quantum dynamics Schollwöck (2011). The inital state of the adiabatic evolution is the band-insulator ground state of the initial lattice for the adiabatic evolution, which is also represented by a matrix-product state in the calculation. To obtain the wave function overlap (see Eq. (S6)) and excitation energy of the final quantum state at the end of the adiabatic evolution, the ground state of the final Hamiltonian is also calculated by the DMRG method. The excitation energy is defined to be energy of the final state of the evolution subtracted by the ground state energy of the final Hamiltonian. The results of wave function overlap and excitation energy are shown in Fig. S5. We see convergence of these results with the bond dimension $D$ increased from $60$ upto $150$. The relative error in excitation energy comparing $D=120$ and $D=150$ is within $2.5\%$ and for the wave function overlap it is within $0.1\%$. We also benchmark our DMRG simulation using the free fermion model with the exact results from the method described in the last section. With the bond dimension $D$ increased to $150$, we find quantitative agreement between the DMRG and the exact results (see Fig. S5(c,d)), which justifies the correctness of the converged DMRG results.