Hybrid-space density matrix renormalization group study of the doped   two-dimensional Hubbard model

G. Ehlers; S. R. White; and R. M. Noack

arXiv:1701.03690·cond-mat.str-el·March 22, 2017

Hybrid-space density matrix renormalization group study of the doped two-dimensional Hubbard model

G. Ehlers, S. R. White, and R. M. Noack

PDF

TL;DR

This paper demonstrates that hybrid real-momentum space DMRG is more efficient than real-space DMRG for the 2D Hubbard model, enabling detailed analysis of charge and magnetic orderings in doped systems.

Contribution

The study introduces a hybrid-space DMRG approach that reduces computational cost and applies it to analyze charge and magnetic orderings in doped 2D Hubbard models.

Findings

01

Hybrid-space DMRG is more efficient than real-space DMRG.

02

Ground state shows stripe charge-density at certain fillings.

03

Charge ordering strength varies with U/t and boundary conditions.

Abstract

The performance of the density matrix renormalization group (DMRG) is strongly influenced by the choice of the local basis of the underlying physical lattice. We demonstrate that, for the two-dimensional Hubbard model, the hybrid real-momentum space formulation of the DMRG is computationally more efficient than the standard real-space formulation. In particular, we show that the computational cost for fixed bond dimension of the hybrid-space DMRG is approximately independent of the width of the lattice, in contrast to the real-space DMRG, for which it is proportional to the width squared. We apply the hybrid-space algorithm to calculate the ground state of the doped two-dimensional Hubbard model on cylinders of width four and six sites; at $n = 0.875$ filling, the ground state exhibits a striped charge-density distribution with a wavelength of eight sites for both $U / t = 4.0$ and $U / t = 8.0$ .…

Tables5

Table 1. Table 1: Average MPO bond dimension for the two-dimensional Hubbard model in the real-space and in the hybrid-space representation for different cylinder widths.

	$L_{y} = 4$	$L_{y} = 6$	$L_{y} = 8$	$L_{y} = 10$
Real space	$18$	$26$	$34$	$42$
Hybrid space	$26.0$	$45.7$	$66.5$	$84.6$

Table 2. Table 2: Scaling of the runtimes of a single Lanczos step, T Lanczos subscript 𝑇 Lanczos {T_{\text{Lanczos}}} , and of a DMRG sweep, T sweep subscript 𝑇 sweep {T_{\text{sweep}}} , and scaling of the memory costs associated with the Krylov space, M Krylov subscript 𝑀 Krylov {M_{\text{Krylov}}} , and the left and right DMRG blocks of the current block-site-site-block configuration, M block subscript 𝑀 block {M_{\text{block}}} , as a function of the cylinder length L x subscript 𝐿 𝑥 L_{x} , the cylinder width L y subscript 𝐿 𝑦 L_{y} , the MPS bond dimension m 𝑚 m , the local lattice dimension d 𝑑 d , and the dimension of the Krylov subspace, K 𝐾 K .

	Real space	Hybrid space
$T_{Lanczos}$	$𝒪 (d^{2} L_{y} m^{3})$	$𝒪 (d^{2} L_{y}^{- 1} m^{3})$
$T_{sweep}$	$𝒪 (d^{2} L_{y}^{2} L_{x} m^{3} K)$	$𝒪 (d^{2} L_{x} m^{3} K)$
$M_{Krylov}$	$𝒪 (d^{2} m^{2} K)$	$𝒪 (d^{2} L_{y}^{- 1} m^{2} K)$
$M_{block}$	$𝒪 (L_{y} m^{2})$	$𝒪 (m^{2})$

Table 3. Table 3: Labels for MPO channels grouped according to the number of included operators as used for the MPO matrix in Eq. ( 45 ).

four		two		one
1	final	6	$𝐜_{↑}^{†}$ $𝐜_{↑}^{}$	12	$𝐜_{↑}^{†}$
three		7	$𝐜_{↑}^{†}$ $𝐜_{↓}^{†}$	13	$𝐜_{↑}^{}$
2	$𝐜_{↑}^{†}$ $𝐜_{↑}^{}$ $𝐜_{↓}^{†}$	8	$𝐜_{↑}^{†}$ $𝐜_{↓}^{}$	14	$𝐜_{↓}^{†}$
3	$𝐜_{↑}^{†}$ $𝐜_{↑}^{}$ $𝐜_{↓}^{}$	9	$𝐜_{↑}^{}$ $𝐜_{↓}^{†}$	15	$𝐜_{↓}^{}$
4	$𝐜_{↑}^{†}$ $𝐜_{↓}^{†}$ $𝐜_{↓}^{}$	10	$𝐜_{↑}^{}$ $𝐜_{↓}^{}$	zero
5	$𝐜_{↑}^{}$ $𝐜_{↓}^{†}$ $𝐜_{↓}^{}$	11	$𝐜_{↓}^{†}$ $𝐜_{↓}^{}$	16	initial

Table 4. Table 4: Zero-truncation-error extrapolated ground-state energies of Hubbard cylinders for U / t = 4.0 𝑈 𝑡 4.0 {U/t=4.0} at n = 0.875 𝑛 0.875 {n=0.875} filling for different system sizes and transverse boundary conditions. For width 6, a pinning field was used to stabilize the different stripe configuration during the initial DMRG sweeps.

$L_{x} \times L_{y}$	boundary conditions	$λ$	$E_{0} / t$
$16 \times 4$	PBC	$8.0$	$- 1.018413 (2)$
$16 \times 4$	ABC	$8.0$	$- 1.001587 (4)$
$32 \times 4$	PBC	$8.0$	$- 1.028614 (6)$
$32 \times 4$	ABC	$8.0$	$- 1.010893 (3)$
$48 \times 4$	PBC	$8.0$	$- 1.032078 (7)$
$48 \times 4$	ABC	$8.0$	$- 1.013997 (5)$
$64 \times 4$	PBC	$8.0$	$- 1.033814 (4)$
$64 \times 4$	ABC	$8.0$	$- 1.015549 (5)$
$\infty \times 4$	PBC	$8.0$	$- 1.03891 (3)$
$\infty \times 4$	ABC	$8.0$	$- 1.020202 (5)$
$16 \times 6$	PBC	$8.0$	$- 1.0097 (5)$
$16 \times 6$	ABC	$8.0$	$- 1.0088 (1)$
$32 \times 6$	PBC	$8.0$	$- 1.0191 (2)$
$32 \times 6$	ABC	$8.0$	$- 1.0184 (1)$
$48 \times 6$	PBC	$8.0$	$- 1.0223 (2)$
$48 \times 6$	ABC	$8.0$	$- 1.0214 (2)$
$64 \times 6$	PBC	$8.0$	$- 1.0240 (3)$
$64 \times 6$	ABC	$8.0$	$- 1.0230 (4)$
$\infty \times 6$	PBC	$8.0$	$- 1.0286 (4)$
$\infty \times 6$	ABC	$8.0$	$- 1.0277 (3)$

Table 5. Table 5: Zero-truncation-error extrapolated ground-state energies of Hubbard cylinders for U / t = 8.0 𝑈 𝑡 8.0 {U/t=8.0} at n = 0.875 𝑛 0.875 {n=0.875} filling for different stripe patterns, system sizes, and transverse boundary conditions. For width 6, a pinning field was used to stabilize the different stripe configuration; the energy contribution of the pinning-field was subtracted afterwards.

$L_{x} \times L_{y}$	boundary conditions	$λ$	$E_{0} / t$
$16 \times 4$	PBC	$8.0$	$- 0.75114 (2)$
$16 \times 4$	ABC	$8.0$	$- 0.74712 (2)$
$32 \times 4$	PBC	$8.0$	$- 0.75841 (2)$
$32 \times 4$	ABC	$8.0$	$- 0.75382 (3)$
$48 \times 4$	PBC	$8.0$	$- 0.76079 (2)$
$48 \times 4$	ABC	$8.0$	$- 0.75604 (4)$
$64 \times 4$	PBC	$8.0$	$- 0.7621 (5)$
$64 \times 4$	ABC	$8.0$	$- 0.75725 (6)$
$\infty \times 4$	PBC	$8.0$	$- 0.7657 (3)$
$\infty \times 4$	ABC	$8.0$	$- 0.76057 (7)$
$16 \times 6$	PBC	$8.0$	$- 0.7481 (2)$
$16 \times 6$	PBC	$5.3$	$- 0.74745 (2)$
$32 \times 6$	PBC	$8.0$	$- 0.7556 (7)$
$32 \times 6$	PBC	$5.3$	$- 0.754702 (3)$
$48 \times 6$	PBC	$8.0$	$- 0.7577 (3)$
$48 \times 6$	PBC	$5.3$	$- 0.75727 (1)$
$64 \times 6$	PBC	$8.0$	$- 0.7591 (2)$
$64 \times 6$	PBC	$5.3$	$- 0.75842 (4)$
$\infty \times 6$	PBC	$8.0$	$- 0.7627 (5)$
$\infty \times 6$	PBC	$5.3$	$- 0.76210 (5)$

Equations49

H = - t ⟨ r, r^{'} ⟩ σ \sum c_{r σ}^{†} c_{r σ}^{†} + U r \sum n_{r ↑}^{†} n_{r ↓}^{†},

H = - t ⟨ r, r^{'} ⟩ σ \sum c_{r σ}^{†} c_{r σ}^{†} + U r \sum n_{r ↑}^{†} n_{r ↓}^{†},

c_{r σ}^{†} = \frac{1}{L _{y}} k_{y} \sum e^{- i k_{y} y} c_{x k_{y} σ}^{†}

c_{r σ}^{†} = \frac{1}{L _{y}} k_{y} \sum e^{- i k_{y} y} c_{x k_{y} σ}^{†}

H =

H =

+ x k_{y} σ \sum ε (k_{y}) n_{x k_{y} σ}^{†}

+ \frac{U}{L _{y}} x k_{y} p_{y} q_{y} \sum c_{x k_{y} + q_{y} ↓}^{†} c_{x p_{y} - q_{y} ↑}^{†} c_{x p_{y} ↑}^{†} c_{x k_{y} ↓}^{†},

H_{n} = P x k_{y} σ \sum sin (ϕ + κ x) n_{x k_{y} σ}^{†}

H_{n} = P x k_{y} σ \sum sin (ϕ + κ x) n_{x k_{y} σ}^{†}

\frac{U}{N} k_{y} \sum A_{x k_{y}}^{†} A_{x k_{y}}^{†}

\frac{U}{N} k_{y} \sum A_{x k_{y}}^{†} A_{x k_{y}}^{†}

A_{x k_{y}}^{†}

A_{x k_{y}}^{†}

A_{x k_{y}}^{†}

S (r, r^{'})

S (r, r^{'})

C (r, r^{'})

D_{y y} (r, r^{'})

S_{z} (x, y)

S_{z} (x, y)

n (x, y)

Δ_{y}^{†} (x, y)

D_{y y} (r, r^{'}) = \frac{1}{2 L _{y}^{2}} k_{y} \sum e^{i k_{y} (y - y^{'})} ⟨ O_{y}^{†} (x, k_{y}) O_{y}^{†} (x^{'}, k_{y})⟩

D_{y y} (r, r^{'}) = \frac{1}{2 L _{y}^{2}} k_{y} \sum e^{i k_{y} (y - y^{'})} ⟨ O_{y}^{†} (x, k_{y}) O_{y}^{†} (x^{'}, k_{y})⟩

O_{y}^{†} (x, k_{y}) = 2 p_{y} \sum cos (p_{y} - k_{y} /2) c_{x p_{y} ↑}^{†} c_{x k_{y} - p_{y} ↓}^{†},

O_{y}^{†} (x, k_{y}) = 2 p_{y} \sum cos (p_{y} - k_{y} /2) c_{x p_{y} ↑}^{†} c_{x k_{y} - p_{y} ↓}^{†},

O_{y}^{†} (x, k_{y}) = 2 p_{y} \sum cos (p_{y} - k_{y} /2) c_{x p_{y} ↑}^{†} c_{x k_{y} - p_{y} ↓}^{†} .

S_{S} (k) = \frac{1}{N} r r^{'} \sum e^{i k (r - r^{'})} S (r, r^{'});

S_{S} (k) = \frac{1}{N} r r^{'} \sum e^{i k (r - r^{'})} S (r, r^{'});

S (l_{x}) = \frac{1}{8} x = (L_{x} - l_{x}) /2 - 4 \sum (L_{x} - l_{x}) /2 + 3 ∣ S (x, x + l_{x}) ∣ .

S (l_{x}) = \frac{1}{8} x = (L_{x} - l_{x}) /2 - 4 \sum (L_{x} - l_{x}) /2 + 3 ∣ S (x, x + l_{x}) ∣ .

c_{↑}^{†} =

c_{↑}^{†} =

c_{↓}^{†} =

1 =

W^{[i]} = 1 t c_{↑}^{†} - t c_{↑}^{†} t c_{↓}^{†} - t c_{↓}^{†} U n_{↑} n_{↓} 00000 c_{↑}^{†} F 00000 c_{↑}^{†} F 00000 c_{↓}^{†} F 00000 c_{↓}^{†} F 000001,

W^{[i]} = 1 t c_{↑}^{†} - t c_{↑}^{†} t c_{↓}^{†} - t c_{↓}^{†} U n_{↑} n_{↓} 00000 c_{↑}^{†} F 00000 c_{↑}^{†} F 00000 c_{↓}^{†} F 00000 c_{↓}^{†} F 000001,

\displaystyle W^{[i]}=\left[{\begin{array}[]{c!{\color{gray}\vrule}cccc!{\color{gray}\vrule}cccc!{\color{gray}\vrule}c}{\bf 1}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr t\,{\bf c}^{\vphantom{{\dagger}}}_{\uparrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ -t\,{\bf c}^{{{\dagger}}}_{\uparrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ t\,{\bf c}^{\vphantom{{\dagger}}}_{\downarrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ -t\,{\bf c}^{{{\dagger}}}_{\downarrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf F}&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&{\bf F}&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&{\bf F}&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&{\bf F}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr\varepsilon_{i}\,{\bf n}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf c}^{{{\dagger}}}_{\uparrow}{\bf F}&{\bf c}^{\vphantom{{\dagger}}}_{\uparrow}{\bf F}&{\bf c}^{{{\dagger}}}_{\downarrow}{\bf F}&{\bf c}^{\vphantom{{\dagger}}}_{\downarrow}{\bf F}\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf 1}\end{array}}\right]\,,}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

\displaystyle W^{[i]}=\left[{\begin{array}[]{c!{\color{gray}\vrule}cccc!{\color{gray}\vrule}cccc!{\color{gray}\vrule}c}{\bf 1}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr t\,{\bf c}^{\vphantom{{\dagger}}}_{\uparrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ -t\,{\bf c}^{{{\dagger}}}_{\uparrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ t\,{\bf c}^{\vphantom{{\dagger}}}_{\downarrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ -t\,{\bf c}^{{{\dagger}}}_{\downarrow}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf F}&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&{\bf F}&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&{\bf F}&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ 0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&{\bf F}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&0\\ \hline\cr\varepsilon_{i}\,{\bf n}\hfil\lx@intercol\color[rgb]{.5,.5,.5&0&0&0&0\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf c}^{{{\dagger}}}_{\uparrow}{\bf F}&{\bf c}^{\vphantom{{\dagger}}}_{\uparrow}{\bf F}&{\bf c}^{{{\dagger}}}_{\downarrow}{\bf F}&{\bf c}^{\vphantom{{\dagger}}}_{\downarrow}{\bf F}\hfil\lx@intercol\color[rgb]{.5,.5,.5&{\bf 1}\end{array}}\right]\,,}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

U / L_{y} k k^{'} p p^{'} \sum c_{k ↑}^{†} c_{k^{'} ↑}^{†} c_{p ↓}^{†} c_{p^{'} ↓}^{†} δ_{k - k^{'} + p - p^{'}},

U / L_{y} k k^{'} p p^{'} \sum c_{k ↑}^{†} c_{k^{'} ↑}^{†} c_{p ↓}^{†} c_{p^{'} ↓}^{†} δ_{k - k^{'} + p - p^{'}},

U / L_{y} k k^{'} p p^{'} \sum c_{k ↑}^{†} c_{k^{'} ↑}^{†} c_{p ↓}^{†} c_{p^{'} ↓}^{†} .

U / L_{y} k k^{'} p p^{'} \sum c_{k ↑}^{†} c_{k^{'} ↑}^{†} c_{p ↓}^{†} c_{p^{'} ↓}^{†} .

W^{[i]} =

W^{[i]} =

\displaystyle\left[\scalebox{0.98}{\mbox{$\displaystyle\begin{array}[]{c|cccc|cccccc|cccc|c}{\bf 1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ \hline\cr U/L_{y}\,{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf F}&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ -U/L_{y}\,{\bf c}^{\dagger}_{\downarrow}&0&{\bf F}&0&0&0&0&0&0&0&0&0&0&0&0&0\\ U/L_{y}\,{\bf c}^{\vphantom{\dagger}}_{\uparrow}&0&0&{\bf F}&0&0&0&0&0&0&0&0&0&0&0&0\\ -U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}&0&0&0&{\bf F}&0&0&0&0&0&0&0&0&0&0&0\\ \hline\cr U/L_{y}\,{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf c}^{\dagger}_{\downarrow}{\bf F}&{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&0&0&{\bf 1}&0&0&0&0&0&0&0&0&0&0\\ -U/L_{y}\,{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&-{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}&0&{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&0&0&{\bf 1}&0&0&0&0&0&0&0&0&0\\ U/L_{y}\,{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}&0&-{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}&-{\bf c}^{\dagger}_{\downarrow}{\bf F}&0&0&0&{\bf 1}&0&0&0&0&0&0&0&0\\ U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf c}^{\dagger}_{\uparrow}{\bf F}&0&0&{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&0&0&0&{\bf 1}&0&0&0&0&0&0&0\\ -U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}&0&{\bf c}^{\dagger}_{\uparrow}{\bf F}&0&-{\bf c}^{\dagger}_{\downarrow}{\bf F}&0&0&0&0&{\bf 1}&0&0&0&0&0&0\\ U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}&0&0&{\bf c}^{\dagger}_{\uparrow}{\bf F}&{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}&0&0&0&0&0&{\bf 1}&0&0&0&0&0\\ \hline\cr{U/L_{y}\,{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}}&{{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf F}}&{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&{{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&0&{\bf c}^{\vphantom{\dagger}}_{\uparrow}&{\bf c}^{\dagger}_{\downarrow}&{\bf c}^{\vphantom{\dagger}}_{\downarrow}&0&0&0&{\bf F}&0&0&0&0\\ {-U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}}&-{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf F}&{-{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&0&{{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&-{\bf c}^{\dagger}_{\uparrow}&0&0&{\bf c}^{\dagger}_{\downarrow}&{\bf c}^{\vphantom{\dagger}}_{\downarrow}&0&0&{\bf F}&0&0&0\\ {U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}}&{{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}}&0&{-{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&-{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&0&-{\bf c}^{\dagger}_{\uparrow}&0&-{\bf c}^{\vphantom{\dagger}}_{\uparrow}&0&{\bf c}^{\vphantom{\dagger}}_{\downarrow}&0&0&{\bf F}&0&0\\ {-U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}}&0&{{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}}&{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf F}&{{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf F}}&0&0&-{\bf c}^{\dagger}_{\uparrow}&0&-{\bf c}^{\vphantom{\dagger}}_{\uparrow}&-{\bf c}^{\dagger}_{\downarrow}&0&0&0&{\bf F}&0\\ \hline\cr U/L_{y}\,{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf F}}&{{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&{{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&{{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}}&{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\uparrow}&{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}&{\bf c}^{\dagger}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\dagger}_{\downarrow}&{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf c}^{\dagger}_{\downarrow}{\bf c}^{\vphantom{\dagger}}_{\downarrow}&{\bf c}^{\dagger}_{\uparrow}{\bf F}&{\bf c}^{\vphantom{\dagger}}_{\uparrow}{\bf F}&{\bf c}^{\dagger}_{\downarrow}{\bf F}&{\bf c}^{\vphantom{\dagger}}_{\downarrow}{\bf F}&{\bf 1}\end{array}$}}\right]\,.

c_{↑}^{†} \to 0000 c_{↑}^{†} 0 00000 c_{↑}^{†} c_{↑}^{†} 00000 0 c_{↑}^{†} 0000 00 c_{↑}^{†} 000 000 c_{↑}^{†} 00 .

c_{↑}^{†} \to 0000 c_{↑}^{†} 0 00000 c_{↑}^{†} c_{↑}^{†} 00000 0 c_{↑}^{†} 0000 00 c_{↑}^{†} 000 000 c_{↑}^{†} 00 .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Hybrid-space density matrix renormalization group study

of the doped two-dimensional Hubbard model

G. Ehlers,1 S. R. White,2 and R. M. Noack1

1 Fachbereich Physik, Philipps-Universität Marburg, 35032 Marburg, Germany

2 Department of Physics and Astronomy, University of California, Irvine, California 92697, USA

Abstract

The performance of the density matrix renormalization group (DMRG) is strongly influenced by the choice of the local basis of the underlying physical lattice. We demonstrate that, for the two-dimensional Hubbard model, the hybrid–real-momentum-space formulation of the DMRG is computationally more efficient than the standard real-space formulation. In particular, we show that the computational cost for fixed bond dimension of the hybrid-space DMRG is approximately independent of the width of the lattice, in contrast to the real-space DMRG, for which it is proportional to the width squared. We apply the hybrid-space algorithm to calculate the ground state of the doped two-dimensional Hubbard model on cylinders of width four and six sites; at ${n=0.875}$ filling, the ground state exhibits a striped charge-density distribution with a wavelength of eight sites for both ${U/t=4.0}$ and ${U/t=8.0}$ . We find that the strength of the charge ordering depends on ${U/t}$ and on the boundary conditions. Furthermore, we investigate the magnetic ordering as well as the decay of the static spin, charge, and pair-field correlation functions.

pacs:

71.10.Fd, 71.27.+a

I Introduction

Although the one-band Hubbard model Hubbard (1963) in two dimensions has long been touted as a leading candidate for explaining the basic phenomenon of high-temperature superconductivity in copper oxide planes Bednorz and Müller (1986), whether the unmodified model provides sufficient features to do this is still controversial Aimi and Imada (2007). In order to clear up this issue, enormous effort is being made to obtain its phase diagram numerically LeBlanc et al. (2015). As the model is doped away from half filling, rich behavior emerges Zheng and Chan (2016), with phases that include antiferromagnetic ordering near half filling, a metallic phase for weak on-site interaction, and a superconducting phase for moderate interaction. Many details of the phase diagram, such as the presence of charge and spin density stripes White and Scalapino (2003); Hager et al. (2005), which have been shown experimentally to play a role in high-temperature superconducting materials Tranquada et al. (1995); Ando et al. (2002), are yet to be unequivocally determined.

One of the most important numerical methods used to address these questions is the density matrix renormalization group (DMRG) White (1992); Schollwöck (2005); Hallberg (2006). In particular, the DMRG can provide unbiased results that are a very useful benchmark and complement for other methods. While the DMRG has been very successful in treating one-dimensional systems and ladders, its application to two-dimensional systems, such as wider cylinders, is much more challenging Stoudenmire and White (2012) because the growth of the entanglement between the DMRG blocks, which is generally proportional to the system width for a short-range Hamiltonian, leads to exponential growth of the computational cost. Furthermore, the DMRG treats two-dimensional systems by mapping the lattice to an intrinsically one-dimensional matrix product state (MPS) Liang and Pang (1994), resulting in longer-range effective interaction even for short-range models.

More recently developed tensor network methods such as multiscale entanglement renormalization (MERA) Vidal (2008); Corboz and Vidal (2009) and projected entangled pair states (PEPS) Verstraete et al. (2008); Corboz et al. (2010) avoid this restriction in principle by adapting the topology of the tensor network to the entanglement structure of the system. However, these methods suffer from the fact that the scaling of their computational costs with the dimension of the Hilbert space treated is much higher than that of the standard DMRG or other MPS-based algorithms.

Another ansatz to broaden the applicability of the DMRG for two-dimensional systems is, instead of changing the topology of the underlying network, to change the local basis of the physical model Krumnow et al. (2016). A change of basis can influence three characteristics that drastically influence the performance of the DMRG: the entanglement between subsystems, the range and structure of the interactions, and the number of quantities that are conserved within the specific representation.

For instance, for translational invariant models, the momentum-space DMRG takes advantage of the conserved momentum quantum number and achieves a significant speedup for a fixed size of the truncated Hilbert space Xiang (1996); Nishimoto et al. (2002). However, the scaling of the size of the Hilbert space needed to maintain a given fixed accuracy is problematic. In particular, the block entropy is zero in the noninteracting limit in the momentum-space representation, as the noninteracting Fermi sea is a product state. Unfortunately, short-range interactions in real space become long-range in momentum space so that all parts of the system become strongly entangled and the block entropy increases rapidly as the interaction is turned on; as a result, the entropy scales with the volume of the system for all nonzero interaction strengths. Thus, the computational cost scales exponentially in the volume of the system, making treating systems of any significant size prohibitively expensive Ehlers et al. (2015).

Recently, it has been shown that, by choosing a mixed basis, one can partially take advantage of the performance benefits of the momentum-space DMRG but also retain the beneficial entanglement scaling of the real-space representation Motruk et al. (2016). Our goal hence is to further explore this approach. In particular, we apply the hybrid–real-momentum-space DMRG to the two-dimensional Hubbard model on a lattice with cylindrical topology.

The remainder of the paper is organized as follows: in Sec. II, we express the Hubbard model in the hybrid-space representation, discuss the structure of the matrix product operator (MPO), and outline how real-space two-point correlation functions can be calculated in hybrid space. In Sec. III, we discuss the computational cost of the hybrid-space DMRG as applied to the two-dimensional Hubbard model. In particular, we analyze the scaling of the CPU-time and memory costs as a function of the cylinder width, and verify the results using realistic numerical calculations. Section IV then describes our study of the ground state of the doped Hubbard model at filling ${n=0.875}$ and interaction strengths ${U/t=4.0}$ and ${U/t=8.0}$ on cylindrical lattices of width $4$ and $6$ . In particular, we address the questions of whether the ground state exhibits stripe structures and whether the pairing correlations are enhanced. Finally, Sec. V contains the conclusion.

II Model and method

II.1 Hubbard model in hybrid space

We investigate the two-dimensional Hubbard model with nearest-neighbor hopping and on-site Coulomb repulsion defined by the Hamiltonian

[TABLE]

where ${\langle{\bf r},{\bf r}^{\prime}\rangle}$ denotes nearest neighbors on a square lattice with lattice sites ${{\bf r}=(x,y)}$ . Here, ${c^{\vphantom{\dagger}}_{{\bf r}\,\sigma}}$ ${(c^{{\dagger}}_{{\bf r}\,\sigma})}$ are creation (annihilation) operators for electrons with spin ${\sigma\in\left\{\uparrow\,,\,\downarrow\right\}}$ , and ${n^{\vphantom{\dagger}}_{{\bf r}\,\sigma}=c^{\dagger}_{{\bf r}\,\sigma}\,c^{\vphantom{\dagger}}_{{\bf r}\,\sigma}}$ is the particle-number operator. We take the lattice geometry to be cylindrical, with cylinder length $L_{x}$ , width (i.e., circumference) $L_{y}$ , a lattice spacing of unity, and periodic or antiperiodic boundary conditions in the transverse direction. Such a geometry is favorable for hybrid–real-momentum-space representation. We Fourier transform in the transverse ( $y$ -) direction, i.e., write

[TABLE]

and correspondingly for ${c{{}^{\vphantom{\dagger}}_{{\bf r}\,\sigma}}}$ . For periodic boundary conditions, the transverse momentum points are given by ${k_{y}=2\pi\,j/L_{y}}$ with integer ${0\leq j<L_{y}}$ , whereas, for antiperiodic boundary conditions, ${k_{y}=2\pi\,(j+\tfrac{1}{2})/L_{y}}$ . The resulting Hamiltonian in hybrid space,

[TABLE]

consists of three terms, a longitudinal hopping term (3a), a transverse hopping term (3b) with dispersion relation ${\varepsilon(k_{y})=-2\,t\cos k_{y}}$ , and a long-range interaction term (3c). Note that the on-site Hubbard interaction becomes long-range in the momentum direction, but remains short-range in the real-space direction.

As will be described in Sec. IV, we find states with striped charge and spin density patterns for systems that are moderately doped away from half filling. In order to stabilize and target states with a particular wavelength of the charge-density stripes, $\lambda$ , we sometimes apply an additional pinning field

[TABLE]

to the doped system. The field couples directly to the local charge density in hybrid space and can be tailored using its amplitude $P$ , wavenumber ${\kappa=2\pi\,/\,\lambda}$ , and phase $\phi$ . Depending on the stability of the target state, we either turn of the pinning field after the initial sweeps of the DMRG, or we keep the field amplitude finite but small throughout the calculation and subtract the contribution of the field to the ground-state energy once the calculation is completed. In order to ascertain that the presence of a small pinning field has no significant effect on the physical results, we have compared calculations with and without a pinning field for cases in which the pinning field can be ramped to zero during the calculation without affecting the final convergence and have found the difference in observables and in energy (with the energy of the pinning field subtracted) to be negligible.

II.2 Hybrid-space matrix product operator

In its modern formulation, the DMRG algorithm is best described within the framework of MPSs and MPOs Östlund and Rommer (1995); Schollwöck (2011). The MPS and MPO store the coefficients of the state ${\left|\Psi\right>}$ and Hamiltonian $H$ as products of matrices associated with each site of the DMRG chain. The MPS and MPO bonds are the contractions of the row and column indices, respectively, of neighboring matrices. The dimensions of the MPS bonds, $m$ , are directly related to the number of states kept in the traditional DMRG, and the dimensions of the MPO bonds refer to the operators stored in the left and right block of the DMRG. The MPO also encodes the rules of how ${H\left|\Psi\right>}$ must be calculated within the DMRG-specific block-site-site-block decomposition of the system and how the DMRG blocks are updated during the sweeping process. Therefore, it is crucial to find the optimal MPO representation for a given model and lattice, i.e., to find the MPO with minimal bond dimension.

Roughly speaking, the dimension of each bond of the MPO depends on the number of terms of the Hamiltonian acting simultaneously on both sites of the bond. If the Hamiltonian is factorizable in an appropriate way, the bond dimension can be minimized by accumulating all interactions between opposing sides of a “composed” operator, so that the same interaction can be expressed using fewer terms Xiang (1996). In this section, we describe how this is done in principle; the details of the construction of the MPO matrices for the hybrid-space Hubbard Hamiltonian are given in Appendix A.

For the two-dimensional Hubbard model in real space, the bond dimension of the optimal MPO is proportional to the system width ${L_{y}}$ . This can be understood as follows:

Figure 1 shows the most common mapping of a two-dimensional square lattice onto the one-dimensional MPS / MPO chain; the one-dimensional path is simply folded over the width of the lattice into two dimensions. Therefore, it is clear that, whenever the chain is cut, the number of nearest-neighbor bonds, and thus the number of hopping term in the two-dimensional Hubbard model in real space, is proportional to the width of the system. Since (almost) all of these terms act on different sites of the lattice, the MPO must include an individual channel for each term; thus its bond dimension is proportional to the system width.

In hybrid space, the situation is more complicated. In particular, the bond dimension of the MPO depends on where we cut the system: if the cut is between two neighboring rings of the cylinder (dashed line A in Fig. 1), only real-space-like hopping terms (3a) are cut, and the bond dimension is again proportional to the system width. If the cut separates sites of the same ring of the cylinder (dashed line B in Fig. 1), the corresponding MPO bond has to carry all the long-range interactions in the term (3c); for each ring this sum runs over three independent momenta. If we were to construct the MPO with individual channels for each term, the resulting bond dimension would be proportional to ${L_{y}^{3}}$ .

Fortunately, for a given bipartition, the sum (3c) can be factorized to reduce the effective number of terms as follows: for fixed $x$ , the part of (3c) for which the annihilation and creation operators are in separated parts of the system can be rewritten as

[TABLE]

with composed operators

[TABLE]

Analogous steps can be carried out for all possible distributions of creation and annihilation operators of both spin species onto the two subsystems, resulting in a formulation in which the total number of terms connecting the two parts of the system is ${\mathcal{O}(L_{y})}$ . Thus, without approximation, we obtain an MPO in which the dimension of each bond is proportional to the system width, just as in real space.

This approach was introduced for the momentum-space DMRG in Ref. Xiang (1996). The technique can be applied to other models as long as the Hamiltonian is factorizable. The hybrid-space Hamiltonian and the optimized MPO for the fermionic Hofstadter model are described in detail in Ref. Motruk et al. (2016). Similar optimizations of MPOs have also been carried out for other systems with long-range interactions, such as the quantum chemical Hamiltonian Keller et al. (2015).

The exact bond dimension for the Hubbard model in hybrid space varies with the position of the cut in the cylinder, whereas in real space it is constant within the bulk of the system.

Assuming a constant MPS bond dimension, a good indicator for the computational cost of the DMRG is the average MPO bond dimension, which is approximately twice as large in the hybrid-space representation as in the real-space case (Table 1). Note that the averaged MPO bond dimension for the hybrid-space MPO may still vary slightly for different orderings of the momentum points of each cylinder ring within the DMRG chain.

II.3 Real-space correlation functions in hybrid space

The equal-time spin, charge, and pair-field correlation functions are defined as

[TABLE]

with

[TABLE]

where ${S_{\rm z}(x,y)}$ measures the local spin, ${n(x,y)}$ is the local charge density, and ${\Delta^{\dagger}_{\rm y}(x,y)}$ is the creation operator for a pair of spin-up and -down particles on sites ${(x,y)}$ and ${(x,y+1)}$ .

Measuring two-point real-space correlation functions in hybrid space raises the same issues as implementing the MPO of the Hamiltonian: applying the Fourier transformation (2) to the correlation functions (7) introduces sums over multiple momenta similar to term (3c). These sums can again be factorized analogously to Eq. (5): the pair-field correlation functions in hybrid space can be written as

[TABLE]

with composed operators

[TABLE]

The spin and charge correlation functions can be treated analogously. Thus, the correlation functions (7) can be measured without changing the scaling of the total computational cost of the algorithm.

III Performance

Here we start by emphasizing that the dimension of the MPS bonds required to obtain a given accuracy scales exponentially with the block entropy; thus, the scaling of the block entropy with system parameters and lattice size has a decisive influence on the performance for fixed accuracy. For two-dimensional systems to which the entropy area law Eisert et al. (2010) applies, the entropy is proportional to system width $L_{y}$ . The resulting exponential scaling of the computational cost with system width is a fundamental limitation of the hybrid-space DMRG as well as of the real-space DMRG and all known MPS-based algorithms.

In this section, we at first neglect the variation in the required bond dimension $m$ on model parameters and lattice size and analyze the scaling of the computational cost of the real-space and hybrid-space DMRG for fixed $m$ . We thus estimate the performance gain of the hybrid-space relative to the real-space representation. We then compare our estimate to measurements of the actual runtime and memory usage for typical calculations and finally come back to the issue of the dependence of the relative accuracy of the hybrid-space and real-space representations on bond dimension $m$ .

III.1 Estimated scaling of the computational cost

The majority of the computational cost in the DMRG algorithm comes from applying the Hamiltonian to a state, ${H|\Psi\rangle}$ ; this is the fundamental step in iterative eigensolvers such as the Lanczos or Davidson algorithms Lanczos (1950); Davidson (1975). Therefore, we examine the scaling of the operations required for this operation with the dimension of the physical lattice sites, $d$ , cylinder width $L_{y}$ , cylinder length $L_{x}$ , and MPS bond dimension $m$ based on the structures of the MPS and MPO used. The computational costs of other operations, e.g., changing the active sites of the DMRG in the sweeping process, benefit from the hybrid-space representation in the same way.

We assume that the bond dimension of the MPO is proportional to $L_{y}$ in both real-space and hybrid-space representations, as is the case for the Hubbard model (see Sec. II.2). In our estimate, we neglect the possibility of exploiting the conservation of spin and charge quantum numbers, as it would have the exact same effect on the computational costs in both representations.

In the real-space representation, one Lanczos step then requires ${\mathcal{O}(d^{2}\,L_{y})}$ multiplications of ${m{\times}m}$ matrices, resulting in ${\mathcal{O}(d^{2}\,L_{y}\,m^{3})}$ operations per Lanczos step and ${\mathcal{O}(d^{2}\,L_{y}^{2}\,L_{x}\,m^{3}\,K)}$ operations per DMRG sweep, with a fixed number of Lanczos steps per DMRG step, $K$ . The corresponding memory costs are ${\mathcal{O}(d^{2}\,m^{2}\,K)}$ for the $K$ -dimensional Krylov space and ${\mathcal{O}(L_{y}\,m^{2})}$ for the left and right DMRG block in the current block-site-site-block configuration.

Next, we analyze the scaling of the costs in the hybrid-space representation under the assumption that the hybrid-space MPS has the same total bond dimension $m$ . Since the Hamiltonian (5) conserves the transverse momentum quantum number, every matrix of the MPS can be written in a block-diagonal form with $L_{y}$ blocks of size ${m^{\prime}{\times}m^{\prime}}$ with ${m^{\prime}\approx m/L_{y}}$ , where every block corresponds to one momentum sector. Therefore, the computational costs reduce to ${\mathcal{O}(d^{2}\,L_{y}^{-1}\,m^{3})}$ per Lanczos step and ${\mathcal{O}(d^{2}\,L_{x}\,m^{3}\,K)}$ per DMRG sweep, and the memory costs become ${\mathcal{O}(d^{2}\,L_{y}^{-1}\,m^{2}\,K)}$ for the Krylov space and ${\mathcal{O}(m^{2})}$ for the DMRG blocks. A side-by-side comparison of all costs is given in Table 2.

Note that the above argument is only valid because the MPS bonds decompose into $L_{y}$ momentum quantum number sectors of approximately equal size ${m^{\prime}}$ . For other conserved quantities such as the total spin or particle number, there typically is a different distribution of the quantum number sectors, with only a few large sectors dominating the computational costs of the algorithm; in these cases, the speedup depends primarily on the size of the largest sectors rather than on the number of sectors.

In conclusion, the computational costs of the hybrid-space DMRG are expected to be independent of $L_{y}$ for fixed $m$ , whereas for the standard real-space DMRG the runtime and the memory consumption scales as ${L_{y}^{2}}$ and $L_{y}$ , respectively (Table 2). Estimating the bulk MPO bond dimension for the Hubbard model to be approximately twice as large in hybrid space as in real space (Table 1), we estimate a total speedup of roughly ${L_{y}^{2}/2}$ .

III.2 Measured performance

In order to investigate the actual performance, we compare the computational costs of real-space and hybrid-space calculations for Hubbard cylinders with length $16$ , widths $4$ , $6$ , and $8$ , and periodic transverse boundary conditions. All calculations were carried out using the same code and using six physical cores on an Intel® Xeon® X5650 CPU. In both representations, we exploit the block-diagonal structure of the MPS and MPO matrices with respect to the conserved charge and total spin quantum numbers; in the hybrid space version, we further decompose the matrices using the transverse momentum quantum number.

Figure 2 shows a comparison of the runtime and memory requirements for different $L_{y}$ as a function of $m$ . As expected, the CPU time per sweep is proportional to $m^{3}$ for large $m$ in both cases, while the peak memory consumption scales with $m^{2}$ [Figs. 2(a) and 2(b)]. The deviation from this limiting behavior for smaller $m$ is due to the quantum number bookkeeping and other overhead in the code and is amplified for the hybrid-space DMRG because of the additional momentum quantum number. In agreement with the predictions above, the speedup of the hybrid-space DMRG over the real-space calculations is larger for wider cylinders [Fig. 2(c)]. Because of the additional overhead, the full speedup of the hybrid-space code is only seen at large $m$ . In the regime in which both methods still provide results within reasonable time in our calculations, the hybrid-space DMRG is approximately $12$ times faster for ${L_{y}=4}$ and up to $20$ and $26$ times faster for ${L_{y}=6}$ and ${L_{y}=8}$ . In terms of the peak memory consumption, we observe a weaker influence of $L_{y}$ : the memory savings for the larger $m$ values treated varies only between a factor of $4$ for ${L_{y}=4}$ and $4.5$ for ${L_{y}=6}$ and ${L_{y}=8}$ [Fig. 2(d)].

In Fig. 3, we compare the CPU time per sweep and the memory consumption of the hybrid-space DMRG for different $L_{y}$ . In agreement with Table 2, the computational costs are almost independent of $L_{y}$ . The minor growth of the runtime that is still present is clearly sublinear in $L_{y}$ . The measured computational costs deviate somewhat from the estimated costs even in the large- $m$ limit, in which the influence of overhead should be insignificant. In particular, the measured absolute runtime of the hybrid-space DMRG is not strictly independent of $L_{y}$ , and thus the observed speedup over real-space DMRG grows more slowly than $L_{y}^{2}$ . This deviation is caused by fact that the MPO dimension is not perfectly proportional to $L_{y}$ for small $L_{y}$ (Table 1).

After having compared the computational costs for fixed $m$ , we now investigate how the change of basis influences the block entropy and the convergence of both methods. Since the MPS bond dimension required to reach a fixed truncation error grows exponentially with the von Neumann entropy of the DMRG blocks, $S(i)$ , even small changes in the entropy can influence the convergence significantly. Figure 4 shows comparisons of $S(i)$ in real and hybrid representations for the half-filled system at ${U/t=4.0}$ , Fig. 4(a), and for the doped system, ${n=0.875}$ , at ${U/t=8.0}$ , Fig. 4(b). In Fig. 4(a), it can be seen that the block entropy in the hybrid-space representation differs only slightly from that in the real-space representation. Note that if we cut the system between neighboring rings of the cylinder (gray lines), the block entropy actually is the same in both cases, as expected. This can also be seen in Fig. 4(b); however, here the entropy in the hybrid representation is perceptibly higher between these points. This illustrates that the nonlocal nature of the interaction within the rings does lead to a moderate increase of the entropy for cuts within the rings in the hybrid representation, especially as ${U/t}$ is increased and the system is doped. Figure 5 illustrates the good agreement of the ground-state energies for fixed $m$ as well as after extrapolation to zero truncation error for both parameter sets in Fig. 4. The slight divergence in the extrapolations for the hybrid versus the real-space calculations, especially evident in Fig. 5(b), shows the limitations of the extrapolation, especially for the real-space results. As can be clearly seen, the maximum MPS bond dimension available is significantly larger for the hybrid-space algorithm, resulting in a more well-controlled extrapolation and higher accuracy in the extrapolated energy. In Fig. 5, we expect that the real-space results would converge towards the hybrid-space results if $m$ could be further increased in the real-space calculation; this is not practically possible here. For narrower systems, the energies converge more rapidly with $m$ , and both methods yield results that agree very well. For example, for the $16{\times}4$ cylinder at $U/t=8.0$ and $n=0.875$ , we obtain the same extrapolated energy, $E_{0}/t=-0.75114(2)$ , for each method.

Even with the significant reductions in computational costs, achieving good convergence for wider cylinders is still very expensive, making an efficient parallelization indispensable. In this respect, the hybrid-space DMRG has an advantageous property: the additional momentum quantum number leads to finer-grained quantum number sectors, which makes for better load balancing. In particular, the momentum quantum numbers yield equal-sized quantum number sectors, resulting in more equal-sized chunks of work. In our implementation, we apply shared-memory parallelization to the combined loops over quantum number sectors and terms of the Hamiltonian, as described for classical DMRG in Ref. Hager et al. (2004). In order to further extend the applicability of the hybrid algorithm, additional steps such as real-space parallelization Stoudenmire and White (2013) or efficient distributed memory parallelization Chan (2004) could be implemented.

IV Results

We study the ground state of the doped Hubbard model on width-4 and width-6 cylinders for ${U/t=4.0}$ and ${U/t=8.0}$ at filling ${n=N^{-1}\sum_{{\bf r}\,\sigma}n_{{\bf r}\,\sigma}=0.875}$ . We choose these two values of $U/t$ for the following reasons: (i) In general, we want to compare the physics and the performance of the method at moderate coupling with that at strong coupling. (ii) It is interesting to investigate the stability of inhomogeneous phases as the interaction strength is changed, in particular, whether a stripe phase remains stable for weaker interaction. (iii) For real-space DMRG methods, the numerical convergence, i.e., the behavior of the block entropy, is generally more favorable at strong coupling. (Real-space DMRG actually becomes exact in the atomic limit rather than the strong-coupling limit, but a local interaction that is strong relative to the hopping generally brings the system closer to that limit.) As was discussed in Sec. III.2, the hybrid-space DMRG has the same block entropy as the real-space method at cuts between cylinder rings and moderate excess entropy for cuts within the rings (see Fig. 4), an excess which becomes smaller for weaker interaction strength. Therefore, the hybrid-space algorithm has slightly better relative convergence at moderate interaction strength and can also handle a larger bond dimension $m$ than the real-space method, so that the moderate-interaction regime is more accessible than with the real-space method.

In the following, we calculate the ground-state energy, the site occupancy, and equal-time spin, charge, and pair-field correlation functions. All results are extrapolated to zero truncation error. The maximum MPS bond dimension was at least ${30\,000}$ for cylinder lengths $16$ and $32$ , ${27\,500}$ for length $48$ , and ${25\,000}$ for length $64$ .

Figure 6 shows the charge-density distribution in the longitudinal direction, ${n(x)=L_{y}^{-1}\sum_{y\,\sigma}n^{\vphantom{\dagger}}_{x\,y\,\sigma}}$ , for ${U/t=4.0}$ and different widths $L_{y}$ and boundary conditions. In all cases, we find a stripe pattern with wavelength 8; for $L_{y}=4$ , each stripe contains $4$ holes [Figs. 6(a) and 6(b)], and for $L_{y}=6$ each stripe contains $6$ holes [Figs. 6(c) and 6(d)]. For larger interaction strength, ${U/t=8.0}$ , Fig. 7, we find two different stripe configurations, wavelength ${\lambda=8.0}$ [Figs. 7(a)–7(c)] and wavelength ${\lambda=5.3}$ (more exactly, ${\lambda=16/3}$ ) [Fig. 7(d)]. For both configurations, we add a pinning field, Eq. (4), with the appropriate wavelength to stabilize the state. In most cases it is sufficient to add the pinning field during the initial sweeps of the calculation and ramp its amplitude $P$ down to zero in subsequent sweeps. However, for width $6$ cylinders at $U/t=8.0$ , the amplitude of the field must be held finite (we take ${P=0.01}$ ) during the entire calculation in order to stabilize the $\lambda=8.0$ stripe pattern. After subtracting the field energy, we find that the ${\lambda=8.0}$ phase is lower in energy and is thus globally stable; it should also be noted that the ${\lambda=5.3}$ phase only occurs for the ${L_{y}=6}$ lattices with periodic boundary conditions.

In order to elucidate the convergence of the $\lambda=5.3$ and $\lambda=8.0$ states for $L_{y}=6$ at $U/t=8.0$ , we examine the discarded-weight extrapolation of both pinning-field-stabilized states in Fig. 8. As can be seen, the $\lambda=8.0$ state has lower energy for any fixed $\Delta\xi$ as well as for $\Delta\xi\to 0.0$ , but the $\lambda=5.3$ state has lower energy for any accessible fixed $m$ . Because of this, the $\lambda=8.0$ state becomes unstable if $P$ is set to zero at any point during the calculation due to the fact that the DMRG algorithm is always driven towards the state that has the lowest energy within the finite reduced subspace of the Hilbert space treated, i.e., the state space with fixed MPS bond dimension $m$ . The ground state within this subspace is not necessarily the ground state for fixed truncation error $\Delta\xi$ or, indeed, for $\Delta\xi\rightarrow 0$ . In general, the DMRG should converge to the true ground state at some level of accuracy (because it becomes exact for $m\rightarrow\infty$ ), but the required $m$ may be inaccessible. Thus, applying a pinning field is necessary here in order to enable the DMRG to track both states as the discarded weight is varied; in practice, we use the fixed pinning field amplitude of $P=0.01$ for the depicted range of $\Delta\xi$ for both states even though it would have been possible to set $P=0$ in the latter part of the $\lambda=5.3$ calculation.

The amplitude of the charge-density modulations, ${\Delta n=1/2\,\{\max_{x}[n(x)]-\min_{x}[n(x)]\}}$ , is plotted as a function of the inverse cylinder length in Fig. 9 for both values of $U/t$ . The trend shows that in all cases the stripes have finite amplitude in the infinite-cylinder length limit. The fact that the stripes are enhanced for ${L_{y}=6}$ indicates a striped ground state in the two-dimensional limit. As the interaction is increased from ${U/t=4.0}$ to ${U/t=8.0}$ , the amplitude of charge density in the ${\lambda=8.0}$ stripes increases.

Typically, a striped charge-density distribution in the doped Hubbard model is accompanied by a striped staggered spin-density distribution with a wavelength that is double that of the charge-density distribution White and Scalapino (2003). Measuring these stripes in hybrid space would require breaking the translation invariance in the transverse direction, which could only be done at the cost of slowing down the algorithm significantly. Instead, we calculate the structure factor of the equal-time spin correlations,

[TABLE]

for finite cylinder length we expand in the particle-in-a-box eigenmodes, as in Ref. Hager et al. (2005). Figures 10 and 11 show ${S_{S}({\bf k})}$ for ${U/t=4.0}$ and ${U/t=8.0}$ , respectively, at ${n=0.875}$ (the same parameter sets as in Figs. 6 and 7). As can be seen in Figs. 10(a)–10(d) for ${U/t=4.0}$ and in Figs. 11(a)–11(c) for ${U/t=8.0}$ , ${S_{S}({\bf k})}$ is peaked at ${{\bf k}\approx(\pm 7/8\,\pi,\pi)}$ in all cases in which period ${\lambda=8.0}$ stripes are present. Thus, the spin correlations are antiferromagnetic with a $\pi$ -phase shift every eighth site, as expected. For the ${\lambda=5.3}$ stripes, Fig. 11(d), the peak shifts to ${{\bf k}\approx(\pm 13/16\,\pi,\pi)}$ which corresponds to antiferromagnetic spin correlations with a $\pi$ -phase shift every $5.3$ sites, also compatible with the stripe structure. As we have seen for the amplitude of the charge-density stripes, the amplitudes of the peaks of the spin structure factor also increase with increasing interaction strength.

In order to test for ${d_{x^{2}-y^{2}}}$ -pairing-induced superconductivity, we have calculated the equal-time pair-field correlation functions and have compared their decay to that of the spin and charge correlations. To compensate for the stripe structure of the ground state, we take averaged absolute values; e.g., for the spin correlation functions we define

[TABLE]

The decay of the correlations with distance along the cylinder is shown in Fig. 12 for width-4 cylinders.

As can be seen by the approximately linear behavior of the envelopes of the correlation functions on the semilogarithmic scale, all three correlation functions decay exponentially at moderate to long distances for both periodic and antiperiodic boundary conditions. For periodic boundary conditions, all three correlation functions decay approximately at the same rate, with the spin correlations having a larger absolute value. For antiperiodic boundary conditions, the spin correlations are dominant for larger distance. Thus, we do not find that pairing correlations are long-range or, indeed, even dominant.

We investigate the effect of interaction strength on the strength of pairing correlations by comparing the pair correlation functions for the ${\lambda=8.0}$ stripes for ${U/t=4.0}$ and ${U/t=8.0}$ in Fig. 13. For periodic boundary conditions, Fig. 13(a), there is a slight suppression of the exponentially decaying correlation function as $U/t$ is increased, whereas for antiperiodic boundary conditions, Fig. 13(b), changing $U/t$ has virtually no effect.

The obvious question to be addressed is why the pair correlations, and, indeed, also the charge and spin correlations, are strong at short length scales, but have no long-range or quasi-long-range (i.e., critical power-law) order. Here we point out that the stripe state breaks the translational symmetry and is locked into the finite lattice, so that the charge order manifests itself in a variation of the local charge density. One then expects the charge correlation with the local order subtracted out to be short-range, i.e., exponentially decaying. A similar argument holds for the spin correlations, which are locked to the charge correlations (but have twice the wavelength and a $\pi$ -phase shift between stripes). The $d_{x^{2}-y^{2}}$ pair correlations do show strong short-range correlations, but are also exponentially decaying, meaning that pairing is not even present in the one-dimensional sense, i.e., that the correlations decay as a power law. To explain this behavior, we make three points: (i) The locked-in charge order could preclude any other than short-range order in all correlation channels. (ii) The pairing correlations are measured perpendicularly to the stripes, i.e., in a direction in which charge transport for static stripes would not be expected. (Note that the cylindrical geometry locks in transverse stripes, and that measurement of pair correlations in the transverse direction over any significant length is not possible due to limitations in the treatable system width.) (iii) For this doping of the system and wavelength of the stripes, the stripes are completely filled with holes and thus insulating, so that quasi-long-range pair correlations would not be expected. In addition, it is fair to point out that, as earlier works on two-leg ladders have shown Noack et al. (1997); Dolfi et al. (2015), convergence of the long-range part of the pairing correlations occurs very slowly as $m$ is increased in the DMRG; an MPS bond dimension of ${30\,000}$ states might still be insufficient to restore algebraic decay over long length scales for width-4 cylinders.

Finally, we obtain accurate estimates of the ground-state energies in the infinite-length-cylinder-limit by carrying out consecutive extrapolations in the discarded weight, followed by the inverse cylinder length.

The extrapolation to infinite cylinder length for ${U/t=4.0}$ and ${U/t=8.0}$ is shown in Fig. 14. Energies for cylinder lengths $16$ , $32$ , $48$ , and $64$ fall almost perfectly onto the linear fits for all curves in both Figs. 14(a) and 14(b). As can be seen, the error bars, which result from the discarded-weight extrapolations, get larger for wider cylinders, but the data points still allow for a well-controlled and accurate extrapolation. Note that in Fig. 14(b) the meta-stable state with a stripe wavelength ${\lambda=5.3}$ is very close in energy to the ${\lambda=8.0}$ state. Even though the effects of the discarded-weight extrapolation, Fig. 8, and the system-length extrapolation can be larger than the energy difference, the energy of the ${\lambda=8.0}$ state is consistently lower than that of the ${\lambda=5.3}$ state, both at finite system length after the discarded-weight extrapolation and after subsequently extrapolating to the infinite-system-length limit. The energies obtained for all parameters are given in Appendix B in Tables 4 and 5.

V Summary and discussion

We have investigated the applicability and usefulness of the DMRG in a hybrid–real-momentum-space formulation for the two-dimensional Hubbard model on cylindrical lattices. We have first compared the computational costs of the real-space and hybrid-space DMRG as a function of the system width in theory and practice, and have shown the hybrid-space variant to be significantly faster. In particular, we have shown that, due; to conservation of the additional momentum quantum number, the computational and memory costs of the hybrid-space DMRG are essentially independent of the width of the cylinder for fixed dimension of the MPS bonds, $m$ . In practice, we have found that the computational cost is almost width-independent and have obtained speedup factors of roughly $12$ , $20$ , and $26$ for cylinder width $4$ , $6$ , and $8$ relative to the real-space DMRG for fixed $m$ . Subsequently, we have demonstrated, that the entropy in the system is not increased significantly in the hybrid representation, giving both formulations comparable convergence with $m$ . Thus, we have shown that the hybrid-space ansatz widens the applicability of the DMRG to Hubbard-like models, especially for larger system widths.

We have then used the hybrid-space DMRG to investigate the static properties of the ground state of the doped two-dimensional Hubbard model at filling ${n=0.875}$ and interaction strengths ${U/t=4.0}$ and ${U/t=8.0}$ ; examining width- $4$ and width- $6$ cylinders, we have found that the system forms a striped charge-density distribution of wavelength ${\lambda=8.0}$ for both values of $U/t$ . The magnetic ordering of the ground state is antiferromagnetic with a modulation of wavelength $16$ . For width 6 cylinders with periodic boundary conditions and interaction strengths ${U/t=8.0}$ , we have also found a metastable ${\lambda=5.3}$ state, which we have shown to be slightly higher in energy for both finite cylinder length and in the infinite-length limit. Furthermore, we have calculated pairing, spin, and charge correlation functions and have found that the spin correlations have the slowest decay with distance, while the pair-field and charge correlations have a sub-dominant decay of comparable strength. The behavior of the correlations is nearly identical for ${U/t=4.0}$ and ${U/t=8.0}$ .

We now compare our results to those of other methods, first for ${U/t=4.0}$ : unsurprisingly, the ground-state energy we have obtained is in excellent agreement with recent real-space DMRG calculations; Leblanc et al. reported a value of ${E_{0}/t=-1.028}$ LeBlanc et al. (2015), for infinite width-6 cylinders averaged over periodic and antiperiodic boundary conditions, while we have obtained ${E_{0}/t=-1.0282(4)}$ . Recent calculations using the density matrix embedding theory (DMET) Knizia and Chan (2012) find a phase diagram in which the system lies just within a superconducting phase for ${n=0.875}$ and ${U/t=4.0}$ . This phase is in very close proximity to both an antiferromagnetic and to a metallic phase; see Fig. 2 in Ref. Zheng and Chan (2016). In particular, the authors find two states with somewhat different magnetic orderings and strengths of $d$ -wave pairing correlations that are energetically very close to one another: an incommensurate antiferromagnet for ${8{\times}2}$ clusters with ${E_{0}/t=-1.0288}$ and a homogeneous antiferromagnetic state for ${4{\times}4}$ clusters with ${E_{0}/t=-1.033(2)}$ . A recent DMET calculation on a ${16{\times}2}$ cluster Zheng and Chan , which is large enough to encompass the $\pi$ -phase shift in spin correlations between stripes, yields a lower energy, ${E_{0}/t=-1.0327(1)}$ , for the incommensurate antiferromagnet. This improved energy for the incommensurate state is within error bounds of that of the homogeneous antiferromagnetic state, leaving the nature of the ground state within the DMET undetermined. The improved energy for the incommensurate antiferromagnetic state and for the homogeneous state agree well with our results; however, we find only a incommensurate striped state, no homogeneous state, in our calculations. We note that we have been able to take larger cluster sizes into account and that an incommensurate striped state is excluded for the ${4{\times}4}$ cluster in the DMET calculations. Therefore, we find it likely that the ground state for this parameter set indeed has an inhomogeneous, i.e., striped, magnetic order in the thermodynamic limit. Earlier DMRG studies of doped width-6 Hubbard cylinders found that the charge-density stripes disappear in the ${L_{x}\rightarrow\infty}$ limit for ${U/t=3.0}$ and ${n\approx 0.905}$ Hager et al. (2005). According to the phase diagram of the DMET calculations, this parameter set falls just on the other side of the phase transition line to a homogeneous metallic state. Thus, these findings are not inconsistent with our results or the DMET results.

The main additional feature of the results for ${U/t=8.0}$ is the presence of a metastable higher energy state with stripe order of wavelength ${\lambda=5.3}$ in addition to the ${\lambda=8.0}$ ground state, found only in the width-6 lattices with periodic boundary conditions. The energy difference between these two states can only be resolved by careful extrapolation in truncation error and system length. We note, however, that these results are consistent with those of other methods, including real-space DMRG, DMET, iPEPS Jordan et al. (2008); Corboz et al. (2014), and constrained-path auxiliary-field quantum Monte Carlo, as discussed extensively in Ref. Zheng et al. (2017). Thus, there is strong evidence that the ${\lambda=8.0}$ stripes are a robust feature of the ground state of the two-dimensional Hubbard model at ${n=0.875}$ for interaction strengths from ${U/t=4.0}$ to ${U/t=8.0}$ .

Acknowledgements.

We thank C.-M. Chung for useful discussions. G.E. and R.M.N. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) through grant no. NO 314/5-2 in Research Unit FOR 1807. S.R.W. was supported by the U.S. National Science Foundation through grant no. DMR-1505406.

Appendix A MPO matrices

The individual matrices of the MPO for the hybrid-space Hamiltonian (3a)–(3c) are relatively large, and their size and structure depend on the width of the cylinder, the mapping between the lattice and the MPO, and the position within the MPO in a nontrivial way. Therefore, it is not possible to give the explicit form of the matrices in a compact way. Instead, we give stepwise instructions on how to construct the explicit form of the MPO.

MPO matrices are best written down as matrices of local operators. For the case of spinful fermions with local basis ${(0,\uparrow,\downarrow,\uparrow\downarrow)}$ , we use

[TABLE]

for the creation, annihilation, identity, and fermionic-sign operators. Here we have incorporated the fermionic sign that occurs when the MPO is applied to an MPS through the ${\bf F}$ operator and the minus signs in ${{\bf c}^{\dagger}_{\downarrow}}$ and ${{\bf c}^{\phantom{\dagger}}_{\downarrow}}$ , which are chosen according to the normal ordering ${({c^{\dagger}_{\uparrow}},{c^{\vphantom{\dagger}}_{\uparrow}},{c^{\dagger}_{\downarrow}},{c^{\vphantom{\dagger}}_{\downarrow}})}$ .

To familiarize the reader with the construction of the fermionic MPO, we start by giving the form of the MPO for the one-dimensional Hubbard model in real space:

[TABLE]

with ${{\bf n}^{\vphantom{{\dagger}}}_{\sigma}={\bf c}^{\vphantom{{\dagger}}}_{\sigma}\,{\bf c}^{{{\dagger}}}_{\sigma}}$ . Here, $W^{[i]}$ is the matrix for the $i$ -th site, and the row and column indices correspond to the virtual MPO bonds. Within each virtual bond, the first and the last channels are used as “initial” and “target” channels of the MPO, and within these channels, identity operators establish connections through the entire MPO from left to right. The on-site repulsion is local in real space and thus its operator directly connects the initial and target channels on each site. For the hopping term, four additional channels are needed to connect the appropriate combinations of creation and annihilation operators on neighboring sites. To complete the MPO, the matrix for the first (last) site must be multiplied with the column (row) unit vector ${\bf e}_{6}^{T}$ ( ${\bf e}_{1}^{\vphantom{T}}$ ).

Using the same basic concepts as in Eq. (14) we can now construct the MPO for the first two parts (3a) and (3b) of the hybrid-space Hamiltonian, i.e., the MPO for the hopping terms. For a width-2 cylinder, the MPO matrices can be written as

[TABLE]

with ${{\bf n}={\bf n}_{\downarrow}+{\bf n}_{\uparrow}}$ . For the mapping between the MPO chain and the two-dimensional lattice given in Fig. 1, an additional four channels are needed to connect the appropriate $c^{\vphantom{\dagger}}_{\sigma}$ and $c^{\dagger}_{\sigma}$ operators between next-nearest-neighbor sites. The transverse hopping is encoded into the dispersion relation and is thus local in hybrid space; it can then be treated just like the on-site repulsion in the real-space MPO. For arbitrary cylinder width, four channels are needed for each momentum point, resulting in a total MPO bond dimension of ${2+4\,L_{y}}$ for the combined terms (3a) and (3b).

The last term (3c) is more complicated to implement, requiring multiple steps. For fixed $x$ , the sum (3c) can alternatively be written as a sum over four momenta,

[TABLE]

where the $\delta$ -function ensures momentum conservation. Note that we have dropped the $x$ index and the second-level $y$ index for compactness. Furthermore, if we neglect momentum conservation for the moment, we obtain a relatively simple sum that contains all possible arrangements of the four operators ${({c^{\dagger}_{\uparrow}},{c^{\vphantom{\dagger}}_{\uparrow}},{c^{\dagger}_{\downarrow}},{c^{\vphantom{\dagger}}_{\downarrow}})}$ within the sites of one ring of the cylinder:

[TABLE]

This sum can be represented as an MPO in a compact form using $16$ channels, which we label in the following way: the initial channel designates “no operator”, the target channel stands for “all four operators”, and the other $14$ channels represent all possible one-, two-, and three-operator combinations (disregarding order and repetition). The complete set of channels is given in Table 3. The lower triangular part of the MPO matrices then contains all entries that logically connect this set of channels in that the matrix element in row $i$ and column $j$ consists of either the operator(s) that must be added to the label of the $i$ -th row to obtain the label of the $j$ -th column or zero if this is not possible. In addition, the matrices must incorporate the appropriate sign that takes into account the fermionic commutation relations with respect to the ordering of the operators in Eq. (28) within the MPO chain. This can be done by assigning a minus sign to all matrix elements that require an odd number of exchanges of operators to reorder the combined labels of the row plus the element itself to match the label of the column. Finally, identity operators and fermionic-sign operators must be placed on the diagonal so that operators connect over longer distances. Following these rules, the resulting MPO matrix for one single ring of the cylinder reads

[TABLE]

Note that fermionic-sign operators must be applied to all elements of the matrix which have an odd number of operators to the right within the MPO, i.e., all elements in columns $2$ to $5$ and $12$ to $15$ . Since the sum (3c) is local in the $x$ -direction, the MPO for the complete cylinder can be constructed by connecting $L_{x}$ MPOs for single rings only through their initial and final channels.

In order to restore momentum conservation, we now split each of the $14$ one-, two-, and three-operator channels into $L_{y}$ separate channels, which correspond to the individual momentum sectors. Only the initial and target channels must have zero momentum and thus do not need to be split. The original matrix (45) must then be transformed in the following way: the corner elements of the matrix remain unaltered, the other elements within the first and last row or column become ${1{\times}L_{y}}$ and ${L_{y}{\times}1}$ sub-matrices, respectively, and all remaining elements become ${L_{y}{\times}L_{y}}$ sub matrices. Within these sub matrices, the local operators must be arranged according to their momentum and the momenta of the column and row channels, i.e., the combined momentum of the row channel and the operator must match the momentum of the column channel (modulo system width). Precisely, if we label the new channels with increasing momenta, an operator with momentum $k$ is placed on the $k$ -th upper and ${(L_{y}-k)}$ -th lower diagonal. For example, a creation operator on a site with momentum ${k=2}$ within a width-6 cylinder becomes a ${6{\times}6}$ sub-matrix

[TABLE]

The resulting MPO is indeed momentum-conserving, and its structure encodes the factorization described in Sec. II.2, with the one-, two-, and three-operator channels corresponding to all possible composed operators to the left and the right side of each bond.

This completes the steps that are necessary to assemble a complete MPO for the hybrid-space Hamiltonian. The MPOs for the different parts of the Hamiltonian can be combined by simply using distinct sets of channels for each part, except for the initial and final channels. In total, one initial channel, one final channel, ${4\,L_{y}}$ channels for the longitudinal hopping (3a), and ${14\,L_{y}}$ channels for the nonlocal Hubbard interaction (3c) are required, resulting in a total virtual bond dimension of ${2+18\,L_{y}}$ .

The final MPO contains many superfluous nonzero elements, which should be eliminated to prevent unnecessary calculations. Since the initial and final channels have zero momentum, many paths through the MPO described above are dead ends that cannot contribute to the final results in any meaningful calculation, i.e., they are either not connected to the initial channel to the left or have no connection to the target channel to the right. One can eliminate all elements that are part of such dead ends and thus significantly reduce the effective average virtual bond dimension of the MPO (see Table 1). Furthermore, it is possible to reuse the one-operator channels of the interaction MPO to implement the longitudinal hopping in a more economical way; however, this makes the structure of the MPO more complicated and results only in minor savings ${(\approx 5\%)}$ of computational time and memory costs.

Appendix B Complete energy tables

Tables 4 and 5 list the ground-state energies for the two-dimensional Hubbard model at ${n=0.875}$ and ${U/t=4.0}$ and ${U/t=8.0}$ , respectively.

Bibliography39

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Hubbard (1963) J. Hubbard, Proc. R. Soc. A 276 , 238 (1963) .
2Bednorz and Müller (1986) J. G. Bednorz and K. A. Müller, Z. Phys. B 64 , 189 (1986) . · doi ↗
3Aimi and Imada (2007) T. Aimi and M. Imada, J. Phy. Soc. Jpn. 76 , 113708 (2007) . · doi ↗
4Le Blanc et al. (2015) J. P. F. Le Blanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng, M. Ferrero, T. M. Henderson, C. A. Jiménez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V. Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White, S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull (Simons Collaboration on the Many-Electron Problem), Phys. Rev. X 5 , 041041 (2015) . · doi ↗
5Zheng and Chan (2016) B.-X. Zheng and G. K.-L. Chan, Phys. Rev. B 93 , 035126 (2016) . · doi ↗
6White and Scalapino (2003) S. R. White and D. J. Scalapino, Phys. Rev. Lett. 91 , 136403 (2003) . · doi ↗
7Hager et al. (2005) G. Hager, G. Wellein, E. Jeckelmann, and H. Fehske, Phys. Rev. B 71 , 075108 (2005) . · doi ↗
8Tranquada et al. (1995) J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature 375 , 561 (1995) . · doi ↗