Quantum state tomography on a plaquette in the 2D Hubbard model

TL;DR

This paper demonstrates that quantum state tomography on a small subsystem within the 2D Hubbard model can be achieved using determinantal quantum Monte Carlo simulations, providing detailed insights into local quantum states.

Contribution

It introduces a method to extract complete quantum state information of a small subsystem from Monte Carlo simulations, including both diagonal and off-diagonal elements of the reduced density matrix.

Findings

Successfully computed the reduced density matrix for a plaquette in the Hubbard model.

Applied symmetry and particle number conservation to analyze eigenvalues.

Provides benchmark data for quantum gas microscope experiments.

Abstract

Motivated by recent quantum gas microscope experiments for fermions in optical lattices, we present proof of principle calculations showing that it is possible to obtain the complete information about the quantum state on a small subsystem from equilibrium determinantal quantum Monte Carlo simulations. Both diagonal (in the occupation number basis) and off-diagonal elements of the reduced density matrix are calculated for a square plaquette, which is embedded in a much larger system of the two-dimensional Hubbard model, both at half filling and in the doped case. The diagonalization of the reduced density matrix is done by exploiting the point group symmetry and particle number conservation, which allows to attach symmetry labels to its eigenvalues. Knowledge of the probabilities of plaquette occupation number configurations is useful for meticulous benchmarking of quantum gas microscope experiments. As the quantum state on the plaquette is exact and self-consistently embedded in an exact, correlated bath, the present approach connects to various cluster approximation techniques.

Tables

Table 1. Table 1: Irreducible representations of the symmetry group D 4 subscript 𝐷 4 D_{4} .

Mulliken symbol	dim.	basis function	symmetry
$A_1$	1	$\left(x^2+y^2\right)\cdot z^2$	$s$
$A_2$	1	$\left(x^2+y^2\right)\cdot z$	$s$
$B_1$	1	$x^2-y^2$	$d_{x^2-y^2}$
$B_2$	1	$xy$	$d_{xy}$
$E$	2	$x,y$	$p_x,p_y$

Table 2. Table 2: Group structure of the Hilbert space for a single square plaquette of the Hubbard model. Shown is the reduction of the subspaces of fixed particle number ( N ↑ , N ↓ ) subscript 𝑁 ↑ subscript 𝑁 ↓ (N_{\uparrow},N_{\downarrow}) into irreducible invariant subspaces of the symmetry group D 4 subscript 𝐷 4 D_{4} . The table is symmetric under exchange of N ↑ subscript 𝑁 ↑ N_{\uparrow} and N ↓ subscript 𝑁 ↓ N_{\downarrow} . Particle number sectors ( N ↑ , N ↓ ) subscript 𝑁 ↑ subscript 𝑁 ↓ (N_{\uparrow},N_{\downarrow}) above half filling have the same group structure as their particle-hole symmetric counterparts ( N s − N ↑ , N s − N ↓ ) subscript 𝑁 𝑠 subscript 𝑁 ↑ subscript 𝑁 𝑠 subscript 𝑁 ↓ (N_{s}-N_{\uparrow},N_{s}-N_{\downarrow}) with N s = 4 subscript 𝑁 𝑠 4 N_{s}=4 .

$N_{\uparrow },N_{\downarrow }$	$\text{Dimension}=\left(\begin{array}{c}{N}_s\\ N_{\uparrow }\end{array}\right)\cdot \left(\begin{array}{c}{N}_s\\ N_{\downarrow }\end{array}\right)$	Irreducible representations
$0,0$	$1$	$A_1$
$1,0$	$4$	$A_1\oplus B_2\oplus E$
$2,0$	$6$	$A_2\oplus B_2\oplus 2E$
$1,1$	$16$	$3A_1\oplus A_2\oplus B_1\oplus 3B_2\oplus 4E$
$3,0$	$4$	$A_2\oplus B_1\oplus E$
$2,1$	$24$	$3A_1\oplus 3A_2\oplus 3B_1\oplus 3B_2\oplus 6E$
$4,0$	$1$	$B_1$
$3,1$	$16$	$A_1\oplus 3A_2\oplus 3B_1\oplus B_2\oplus 4E$
$2,2$	$36$	$6A_1\oplus 4A_2\oplus 6B_1\oplus 4B_2\oplus 8E$

Table 3. Table 3: Group structure of the spin-balanced subspaces of the single-plaquette Hubbard model. The decomposition of the subspaces is done with respect to the irreducible invariant subspaces of D 4 ​ h subscript 𝐷 4 ℎ D_{4h} . The subscript g 𝑔 g ( u 𝑢 u ) denotes a wave function that is even (odd) under flipping all spins on the plaquette.

$N_{\uparrow },N_{\downarrow }$	Irreducible representations
$1,1$	$3A_{1,g}\oplus B_{1,g}\oplus 2B_{2,g}\oplus 2E_g\oplus A_{2,u}\oplus B_{2,u}\oplus 2E_u$
$2,2$	$5A_{1,g}\oplus A_{2,g}\oplus 4B_{1,g}\oplus 3B_{2,g}\oplus 4E_g\oplus A_{1,u}\oplus 3A_{2,u}\oplus 2B_{1,u}\oplus B_{2,u}\oplus 4E_u$

Equations

$$\mathcal{H}=-$$

$$+$$

$$Z=\text{Tr}\left(e^{-\beta \mathcal{H}}\right)$$

$$Z$$

$$=\sum _{\{\mathbf{s}\}}{w}_{\{\mathbf{s}\}}^{\uparrow }{w}_{\{\mathbf{s}\}}^{\downarrow }.$$

$${\left[G_{\{\mathbf{s}\}}^{\sigma }(\tau =l\Delta \tau )\right]}_{i,j}\equiv ⟨c_{i,\sigma }{c}_{j,\sigma }^{†}{⟩}_{\{\mathbf{s}\}}$$

$$={\left[{\left(1+B_l^{\sigma }{B}_{l-1}^{\sigma }\cdots B_1^{\sigma }{B}_{N_{\tau }}^{\sigma }\cdots B_{l+1}^{\sigma }\right)}^{-1}\right]}_{i,j},$$

$$⟨\alpha |\rho |\alpha ⟩$$

$$⟨\beta |\rho |\alpha ⟩$$

$$|\alpha ⟩=|{\alpha }_{\uparrow }⟩\otimes |{\alpha }_{\downarrow }⟩=|n_{1,\uparrow },n_{1,\downarrow };n_{2,\uparrow },n_{2,\downarrow };\dots ;n_{N_s,\uparrow },n_{N_s,\downarrow }⟩$$

$${⟨\left(c_{i_1}{c}_{j_1}^{†}\right)\left(c_{i_2}{c}_{j_2}^{†}\right)\cdots \left(c_{i_n}{c}_{j_n}^{†}\right)⟩}_0=\det \left(G_{i_{\alpha }{j}_{\beta }}^{(0)}\right)$$

$$\begin{array}{r}\left(⟨{\beta }_{\uparrow }|\otimes ⟨{\beta }_{\downarrow }|\right){\rho }_{A,\{\mathbf{s}\}}^{\uparrow }\otimes {\rho }_{A,\{\mathbf{s}\}}^{\downarrow }\left(|{\alpha }_{\uparrow }⟩\otimes |{\alpha }_{\downarrow }⟩\right)\\ =⟨{\alpha }_{\uparrow }|{\rho }_{A,\{\mathbf{s}\}}^{\uparrow }|{\beta }_{\uparrow }⟩⟨{\alpha }_{\downarrow }|{\rho }_{A,\{\mathbf{s}\}}^{\downarrow }|{\beta }_{\downarrow }⟩,\end{array}$$

$${\hat{\Pi }}_{\alpha }=\prod _{i\in N_s}\left[{\hat{n}}_i{n}_i+(1-{\hat{n}}_i)(1-n_i)\right]$$

$$=\prod _{i\text{ occupied}}\widehat{n_i}\prod _{j\text{ unoccupied}}(1-{\hat{n}}_j).$$

$${\hat{\Xi }}_{\alpha \to \beta }={\hat{T}}_{\alpha \to \beta }{\hat{\Pi }}_{\alpha },$$

$${\hat{T}}_{\alpha \to \beta }=(-1{)}^p\left(\prod _{\begin{array}{c}c=N^{+}\\ i_1<i_2<\cdots <i_{N^{+}}\end{array}}^1{c}_{i_c}^{†}\right)\left(\prod _{\begin{array}{c}a=N^{-}\\ j_1<j_2<\cdots <j_{N^{-}}\end{array}}^1{c}_{j_a}\right).$$

$$(-1{)}^p=\prod _{i_c=1}^{N^{+}}{e}^{i\pi {\sum }_{i_c<l<N_s}{n}_l^{(\gamma )}}\prod _{j_a=1}^{N^{-}}{e}^{i\pi {\sum }_{j_a<k<N_s}{n}_k^{(\alpha )}}$$

$$(-1{)}^{p^{\prime }}=(-1{)}^{{\sum }_{i=1}^{N^{+}}i}=(-1{)}^{\frac{N^{+}}{2}(N^{+}+1)}.$$

$$\rho =\frac{1}{Z}{e}^{-\beta (H-\mu N)}\text{with}Z=\text{Tr}\left(e^{-\beta H}\right)$$

$$U^{†}\rho U=\rho .$$

$$D_4=\{\mathcal{E},C_{2z},C_{4z},C_{4z}^{-1},C_{2x},C_{2y},C_{2xy},C_{2x\bar{y}}\}$$

$$R=R^{(A)}{R}^{(B)}\text{with}{R}^{(A)}{\mathcal{L}}_A={\mathcal{L}}_A\text{and}{R}^{(B)}{\mathcal{L}}_B={\mathcal{L}}_B,$$

$$P_R=P_R^{(A)}\otimes P_R^{(B)}$$

$$P_R^{†}\rho P_R=\rho$$

$${\text{Tr}}_B\left({P_{R^{\prime }}^{(A)}}^{†}\otimes {P_{R^{\prime }}^{(B)}}^{†}\rho P_{R^{\prime }}^{(A)}\otimes P_{R^{\prime }}^{(B)}\right)$$

$$\Rightarrow {P_{R^{\prime }}^{(A)}}^{†}{\text{Tr}}_B\left({P_{R^{\prime }}^{(B)}}^{†}\rho P_{R^{\prime }}^{(B)}\right)P_{R^{\prime }}^{(A)}$$

$$\Rightarrow {P_{R^{\prime }}^{(A)}}^{†}{\rho }_A{P}_{R^{\prime }}^{(A)}$$

$${\rho }_A^{(\text{representation})}=S^{†}{\rho }_A^{(\text{occupation})}S$$

$$⟨2|{\Delta }_d|4⟩\ne 0$$

$$c_{i,\downarrow }⟶(-1{)}^{\mathbf{i}}{c}_{i,\downarrow }^{†},c_{i,\uparrow }⟶c_{i,\uparrow },$$

$$(2,2)B_{1,g}\ni |{\psi }_1⟩\sim \frac{1}{\sqrt{2}}(|\downarrow \uparrow \uparrow \downarrow ⟩+|\uparrow \downarrow \downarrow \uparrow ⟩)+\dots$$

$$(2,2)A_{2,u}\ni |{\psi }_2⟩\sim \frac{1}{\sqrt{2}}(|\downarrow \uparrow \uparrow \downarrow ⟩-|\uparrow \downarrow \downarrow \uparrow ⟩)+\mathcal{O}(t/U).$$

$$\begin{array}{r}(1,1)A_{1,g}\ni |{\psi }_3⟩\sim \frac{1}{2}(|\uparrow hh\downarrow ⟩+|h\uparrow \downarrow h⟩\\ +|h\downarrow \uparrow h⟩+|\downarrow hh\uparrow ⟩)+\dots ,\end{array}$$

$$\begin{array}{r}(1,3)A_2\ni |{\psi }_4⟩\sim \frac{1}{2}(|\uparrow \downarrow \downarrow \downarrow ⟩+|\downarrow \uparrow \downarrow \downarrow ⟩\\ -|\downarrow \downarrow \uparrow \downarrow ⟩-|\downarrow \downarrow \downarrow \uparrow ⟩)+\mathcal{O}(t/U)\end{array}$$

Full text

Knowledge of the probabilities of plaquette occupation number configurations is useful for meticulous benchmarking of quantum gas microscope experiments. As the quantum state on the plaquette is exact and self-consistently embedded in an exact, correlated bath, the present approach connects to various cluster approximation techniques.

I Introduction

Quantum state tomography refers to the task of reconstructing the full quantum state of a system from measurements, which, by definition, is a task that scales exponentially in the system size. Tomography has been used to characterize small quantum systems consisting of a few qubits, such as trapped ion chains Häffner et al. (2005), molecules in nuclear magnetic resonance experiments Xin et al. (2017), superconducting circuits Baur et al. (2012), and photonic systems Schwemmer et al. (2014); Gao et al. (2018). More efficient quantum tomography techniques based on compressed sensing Gross et al. (2010) and matrix product states Cramer et al. (2010) have been proposed, facilitating the reconstruction based on incomplete data Schwemmer et al. (2014); Riofrío et al. (2017); Steffens et al. (2015).

For cold atoms in optical lattices emulating strongly correlated solid state models such a global characterization of the state is neither feasible nor meaningful. Still, the local state on a subsystem of sites, encoded by its reduced density matrix $\rho_{A}=\text{Tr}_{B=\bar{A}}\left(\rho\right)$, which is given by the partial trace of the global density matrix $\rho=e^{-\beta\hat{H}}/Z$ over the complement system $B$, can provide valuable information if subjected to a tomographic measurement. Here, we present proof of principle determinantal quantum Monte Carlo (DQMC) calculations for a square plaquette of the Fermi-Hubbard model which is embedded in a much larger system.

Cold atoms experiments with quantum gas microscopes have reached a regime where local antiferromangetic correlations in the two-dimensional Fermi-Hubbard model and the effect of doping away from half filling can be explored Mazurenko et al. (2017); Cheuk et al. (2016); Mitra et al. (2017); Koepsell et al. (2018). In particular, signatures of polarons Koepsell et al. (2018); Chiu et al. (2018) have been reported.

The finite-temperature phase diagram of the 2D Fermi Hubbard model for intermediate interaction strength Šimkovic et al. (2017) is very challenging due to the fermionic sign problem, which is aggravated exponentially by decreasing temperature. The temperatures currently realized in fermionic quantum gas microscopes experiments are still within reach of numerically exact DQMC simulations, and detailed comparison of full particle number distribution functions Humeniuk and Büchler (2017) has been made.

There are numerous ramifications that motivate computing the equilibrium probabilities of individual microstates, i.e. the diagonal elements of $\rho_{A}$. By superimposing a specific experimental measurement protocol, the effect of parity projection Parsons et al. (2016), i.e. the inability to dinstinguish doubly occupied sites and holes in the imaging process, on the experimentally observed particle configurations can be assessed.

Remarkably, a scheme for measuring the second moment of the density matrix, its purity $\text{Tr}\left(\rho_{A}^{2}\right)$, without exponential effort has been demonstrated in bosonic cold atoms experiments Daley et al. (2012); Islam et al. (2015) and generalized to fermionic systems Pichler et al. (2013). Measurement protocols for accessing the full entanglement spectrum in cold atoms systems have also been proposed Pichler et al. (2016). In view of this and recent progress in using machine learning with neural networks to reconstruct a full quantum state of bosonic systems 111The sign structure of fermionic wave functions makes the application of machine learning approaches to fermionic systems more involved. from a limited number of experimental measurements Torlai et al. (2018, 2019), numerical access to all (diagonal and off-diagonal) elements of $\rho_{A}$ is of potential interest. Alternative numerical approaches relying on the replica representation of Rényi entropies Chung et al. (2014) $\text{Tr}\left[\rho_{A}^{n}\right]$ require complicated modifications in the topology of the simulation cell, whereas the brute force numerical scheme described here does not affect the core of the DQMC algorithm and all diagonal and off-diagonal elements of the reduced density matrix, as well as its eigenvalues, the entanglement spectrum, can be obtained. However, computational and memory resources that grow exponentially with the subsystem size limit the latter to maximally $N_{s}=9$ sites 222The computational complexity of DQMC for simulating the total system of $N$ sites is $\sim\beta N^{3}$ and can be completely decoupled from the costly “exact diagonalization” inside each HS sample if the single-particle Green’s functions are saved on disk for every HS configuration (or after a number of Monte Carlo steps proportional to the autocorrelation time). This requires several hundred GB of hard disk memory per parameter set $(\beta,U)$..

From a methodological point of view, there are connections with various types of numerical cluster approaches. Eqs. (6) and (7) below give the exact state on a cluster that is self-consistently embedded in a correlated bath and can be used to compare with computationally less expensive methods that solve the cluster system exactly, but treat the bath degrees of freedom only approximately. In Refs. Udagawa and Motome (2010, 2015), using an auxiliary-field quantum Monte Carlo solver in the context of cellular dynamical mean-field theory, the entanglement spectrum was computed for a triangular plaquette of the kagome Hubbard model, revealing an emergent composite degree of freedom due to geometric frustration. The method used there is very similar in spirit to the one presented in this work.

Finally, a strong motivation for studying the detailed structure of a local quantum state comes from the phenomenology of the high-temperature phase of the repulsive Hubbard model (or more generally of high-$T_{c}$ superconductors) where for temperatures $\beta\leq 4-5$, where DQMC simulations are still possible due to a mild sign problem Iglovikov et al. (2015), a pseudogap develops in the single-particle spectral function Preuss et al. (1995). The pseudogap in the attractive Hubbard model is well understood in terms of local bound pairs of fermions without long-range phase coherence, the gap being associated with the binding energy of the pair. A natural question is whether similar preformed objects are responsible for the pseudogap observed in the repulsive Hubbard model Preuss et al. (1995) (or more generally in the normal state of high-$T_{c}$ superconductors) and which signatures of correlated phases are contained in the local density matrix Huber et al. (2019).

The structure of this paper is as follows. In Sect. II, we describe the algorithm for projecting the reduced density matrix $\rho_{A}$ on a subsystem $A$ from the global density matrix which is sampled in the DQMC procedure. The symmetries of the reduced density matrix on a square plaquette and the transformation to the irreducible representation basis is discussed in Sect. III. In Sect. IV the issue of error bars is addressed. Finally, Sect. V contains the results for the tomographic reconstruction of $\rho_{A}$ on a square plaquette of the Hubbard model: At half filling, we present the diagonal and off-diagonal elements of $\rho_{A}$ as a function of Hubbard interaction $U/t$ for both high and low temperatures. Away from half-filling, where the computational cost of DQMC simulations is affected by the sign problem, we show the doping dependence of $\rho_{A}$ for $U/t=6$ and a relatively high temperature of $T=0.25t$, which corresponds to the lowest temperature achieved so far in fermionic cold atoms experiments Mazurenko et al. (2017). Sect. VI concludes with an outlook.

II Born’s rule for many-body states

The Hamiltonian studied in the following is that of the single-band Hubbard model

[TABLE]

where $c_{i,\sigma}^{\dagger}$ creates a fermion with spin $\sigma$ at site $i$, and $n_{i,\sigma}=c_{i,\sigma}^{\dagger}c_{i,\sigma}$. Here, $t$ is the hopping matrix element between nearest neighbours $\langle i,j\rangle$, $\mu$ is the chemical potential, and $U$ is the on-site repulsive ($U>0$) or attractive ($U<0$) interaction. The partition sum

[TABLE]

is sampled with the determinantal quantum Monte Carlo method Blankenbecler et al. (1981); Loh Jr. and Gubernatis (1992); Assaad . We briefly sketch the essential conceptual steps in the derivation of this procedure, referring to the exhaustive literature (see e.g. Refs. Loh Jr. and Gubernatis (1992); Assaad ) for more details. After discretizing inverse temperature $\beta=N_{\tau}\Delta\tau$ into $N_{\tau}$ imaginary time slices and separating the single-body kinetic term in the Hamiltonian from the two-body interaction term via a Trotter-Suzuki decomposition, a Hubbard-Stratonovich (HS) transformation is applied to the interaction term converting it into a single-particle term which is coupled to a fluctuating space- and imaginary time-dependent potential given by auxiliary field variables. Thanks to the HS transformation the partition sum contains only exponentials of bilinear (i.e. free) fermionic operators. Then, the free fermions can be integrated out for each auxiliary field configuration using the well-known formula for the grand-canonical fermionic trace, which results in

[TABLE]

Here, $B_{l}^{\sigma}\equiv e^{-\Delta\tau V_{l}^{\sigma}(\{{\bf s}_{l}\})}e^{-\Delta\tau K}$ is an $N\times N$ matrix representation of the single-particle propagator at time slice $l$ after HS transformation, with $V_{l}^{\sigma}(\{{\bf s}_{l}\})$ denoting the potential term and $K$ the kinetic term of the single-particle Hamiltonian after HS transformation Assaad , while $\{{\bf s}\}\equiv\{s_{i,l}\}_{i=1,\ldots,N;l=1,\ldots,N_{\tau}}$ is the space-time configuration of auxiliary-field variables, which is sampled with a Monte Carlo technique.

It can be shown (see e.g. Ref. Assaad ) that the weight $w^{\sigma}_{\{{\bf s}\}}$ is given by the inverse determinant of the equal-time single-particle Green’s function

[TABLE]

which constitutes the central quantity of the DQMC algorithm. Furthermore, for the Hubbard model the HS transformation can be chosen such that the weight of an auxiliary-field configuration $\{{\bf s}\}$ factorizes between spin species Assaad .

Thus, loosely speaking, the DQMC framework consists in computing a large sum over free fermions systems in varying external potentials: It suffices to compute any quantity (in any single-particle basis) for free fermions and average over Monte Carlo samples. This is a huge conceptual simplification compared to path integral methods. The generality of the free fermion decomposition Grover (2013) allows to carry the measurement part of the DQMC algorithm to its extreme by calculating the full quantum state on a small subsystem $A$, that is all elements of the reduced density matrix $\rho_{A}$. This amounts to performing exact diagonalization inside the measurement part of the Monte Carlo procedure.

The diagonal and off-diagonal elements of the reduced density matrix in the Fock basis can be written as

[TABLE]

Here, $\hat{\Pi}_{\alpha}=|\alpha\rangle\langle\alpha|$ are projectors onto individual Fock states

[TABLE]

on the substem $A$ with $N_{s}$ sites 333Translational invariance can be used to accumulate additional statistics by displacing the subsystem repeatedly. However, for the calculations presented in the following for a single plaquette, the location of the plaquette was fixed. and $\hat{\Xi}_{\alpha\rightarrow\beta}=|\beta\rangle\langle\alpha|$ is a transition operator between two Fock states $|\alpha\rangle$ and $|\beta\rangle$. $\rho_{\{\bf s\}}$ is the global density matrix of the free fermion system with auxiliary field configuration $\{{\bf s}\}$, which is sampled via Monte Carlo.

In a state of a non-interacting Fermi system, Wick’s theorem applied to a product of $n$ pairings of fermionic operators results in the determinant formula

[TABLE]

with $\alpha,\beta=1,\ldots,n$, where the equal-time Green’s function of the non-interacting Fermi system is $G_{i_{\alpha}j_{\beta}}^{(0)}(\tau=0)=\langle c_{i_{\alpha}}(\tau=0)c_{j_{\beta}}^{\dagger}(0)\rangle_{0}$. This formula is the basis for evaluating the expectation values of the projectors and transition operators in Eqs. (6) and (7) and thus extracting the elements of the reduced density matrix from an equilibrium state, which in the DQMC framework is encoded in a sum over free fermion systems parametrized by auxiliary field configurations $\{{\bf s}\}$.

Since for one auxiliary field configuration $\{{\bf s}\}$ the reduced density matrix factorizes between spin species,

[TABLE]

at most $2\times 2^{N_{s}}\times 2^{N_{s}}$ matrix elements need to be computed to express all $4^{N_{s}}\times 4^{N_{s}}$ elements of $\rho_{A,\{{\bf s}\}}$. (In this crude estimate we have disregarded the block diagonal structure of $\rho_{A,\{{\bf s}\}}^{\sigma}$ with respect to particle number $N^{\sigma}_{A}$ which reduces the size of the largest particle number block for one spin species to $\begin{pmatrix}N_{s}\\ \lfloor N_{s}/2\rfloor\end{pmatrix}\cdot\begin{pmatrix}N_{s}\\ \lfloor N_{s}/2\rfloor\end{pmatrix}$ with $\lfloor x\rfloor$ denoting the largest integer that is smaller than $x$.) Therefore, the limiting factor is the memory requirement for storing all elements of $\rho_{A,\{{\bf s}\}}$ for Monte Carlo averaging, rather than the computation of individual elements. The presence of point group symmetry operations which leave subsystem $A$ invariant or spin inversion symmetry leads to an additional block diagonal structure of $\rho_{A}$ (see Sect. III.1).

In the following, we discuss the computation for a single spin component, thereby dropping all spin indices in the notation. We use hats to distinguish the number operator $\hat{n}_{i}=c_{i}^{\dagger}c_{i}$ from the occupation number $n_{i}$ and write

[TABLE]

The transition operator $\hat{\Xi}_{\alpha\rightarrow\beta}$ can be written as

[TABLE]

where $\hat{\Pi}_{\alpha}$ projects onto the Fock state $|\alpha\rangle$, which is then converted into $|\beta\rangle$ by a combination of appropriately chosen creation and annihilation operators

[TABLE]

The sequences of site indices $\mathcal{I}^{+}=\{i_{1},i_{2},\ldots,i_{N^{+}}\}$ and $\mathcal{I}^{-}=\{j_{1},j_{2},\ldots,j_{N^{-}}\}$, ordered according to the chosen fermion ordering, denote the lattice sites where the $N^{+}$ creation and $N^{-}$ annihilation operators must act to convert $|\alpha\rangle$ into $|\beta\rangle$. Since $\rho_{A}$ is block diagonal with respect to the total particle number $N_{A,\sigma}=\sum_{i\in A}n_{i,\sigma}$ for each spin species $\sigma$ (see Sect. III.1), there must be as many creation as annihilation operators and $N^{+}=N^{-}$. The fermionic phase

[TABLE]

with ${\bf n}^{(\alpha)}$ the vector of occupation numbers for state $|\alpha\rangle$ and ${\bf n}^{(\gamma)}$ for state $|\gamma\rangle\equiv\prod_{a=N^{-}}^{1}\,c_{j_{a}}|\alpha\rangle$ ensures that $\hat{T}_{\alpha\rightarrow\beta}|\alpha\rangle=|\beta\rangle$. Given that $N^{+}=N^{-}$, one can bring the operator product in Eq. (13) into the paired form as it appears on the left-hand side of Eq. (9) by commuting fermionic operators. From this an additional phase factor arises:

[TABLE]

With regard to a practical implementation for the evaluation of the expectation values in the right-hand side of Eqs. (6) and (7) for one particular auxiliary field configuration $\{{\bf s}\}$, a few remarks are in order: The fact that an occupied site (“occ.” in the code listing in Appendix D) is represented by a projector of the form $\hat{n}_{i}=(1-c_{i}c^{\dagger}_{i})$ (for unoccupied sites $1-\hat{n}_{i}=c_{i}c^{\dagger}_{i}$) means that the expectation value of the total projector onto a Fock state Eq. (11), when multiplied out, is a sum of terms which can be written as a binary tree for the occupied sites where each leaf is of the form of Eq. (9). The branches of the binary tree need to be summed over to obtain the projector Eq. (11). In combination with $\hat{T}_{\alpha\rightarrow\beta}$, the projector $\hat{\Pi}_{\alpha}$ needs only be realized on sites that are unaffected (“ua” in the code listing in Appendix D) by the hopping operators in $\hat{T}_{\alpha\rightarrow\beta}$ since the hopping operators already guarantee that occupation number states $\alpha^{\prime}\neq\alpha$ are eliminated by the action of $\hat{T}_{\alpha\rightarrow\beta}$.

The algorithm for computing $\langle\beta|\rho_{A}|\alpha\rangle_{\{{\bf s}\}}$ between the occupation number states $|\alpha\rangle$ and $|\beta\rangle$ for a single Hubbard-Stratonovich configuration $\{{\bf s}\}$ is summarized in the pseudocode listing in Appendix  D, where the main task consists in collecting the correct row and column indices for the submatrices that enter the determinant formula Eq. (9). The final result for a matrix element of $\rho_{A}$ in the interacting system is obtained by summing over all Hubbard-Stratonovich configurations $\{{\bf s}\}$. As an illustration, the Monte Carlo timeseries of the diagonal element of $\rho_{A}$ corresponding to one of the two Néel states on a plaquette is displayed in Fig. 1. There are rare outliers that exceed to maximum probability of $1$. The histogram of probabilities on the right is slightly bimodal but smooth; the mean value is around $P_{\text{N\'{e}el}}\approx 0.12$. For a projective measurement one would expect a binary distribution of probabilities with only the probabilities [math] or $1$ appearing.

III Symmetries of the reduced density matrix

The form of the global density matrix

[TABLE]

implies that a symmetry operation represented by a unitary operator $U$ which obeys the commutation relation $[H-\mu N,U]=0$, is trivially also a symmetry of the global density matrix:

[TABLE]

In the following it is discussed how symmetries of the Hamiltonian affect the block diagonal structure of the reduced density matrix $\rho_{A}$ of a subsystem $A$. Details on the exact diagonalization of the Hubbard model by means of symmetries are discussed in Refs. Fano et al. (1992); Noce and Cuoco (1996); an analytical diagonalization exploiting all symmetries was carried out in Ref. Schumann (2002) for a single plaquette of the Hubbard model and in Ref. Kuns et al. (2011) for a plaquette of the $t-J$ model.

III.1 Particle number conservation

The eigenstates of the reduced density matrix retain their good quantum numbers when the corresponding operator of the total system is a direct sum of the operators of its subsystems Chen and Moukouri (1996). This is the case for the particle number $\hat{N}_{\sigma}=\hat{N}_{A,\sigma}+\hat{N}_{B,\sigma}$ (or the $z$ component of the total spin $\hat{M}=\hat{N}_{\uparrow}-\hat{N}_{\downarrow}$) and consequently $\rho_{A}$ is block diagonal with respect to the quantum numbers $N_{A,\sigma}=\sum_{{\bf i}\in A}n_{{\bf i},\sigma}$. For brevity, we denote particle number sectors on the subsystem $A$ in the following by $(N_{\uparrow},N_{\downarrow})$, where $N_{\sigma}\equiv N_{A,\sigma}$.

III.2 Lattice symmetry: point group $D_{4}$

We focus in the following on the point group symmetry of the square lattice, the non-Abelian dihedral group $D_{4}$ with $h=8$ group elements

[TABLE]

comprising the identity $\mathcal{E}$ and (assuming that the square is lying in the $x-y$ plane) rotations by $\pi$ around the $x$, $y$ and $z$ axes, $C_{2x},C_{2y}$ and $C_{2z}$, rotations by $\pi$ around the diagonal lines $x=y$ and $x=-y$, $C_{2xy}$ and $C_{2x\bar{y}}$, and clockwise and counterclockwise rotations by $\pi/2$ around the $z$ axis, $C_{4z}$ and $C_{4z}^{-1}$.

The group $D_{4}$ has five irreducible representations, four one-dimensional representations with the Mulliken symbols $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$, and one two-dimensional representation $E$. For later reference, they are listed in Tab. 1 together with their symmetries. (see Tinkham (1964) for the character table and the irreducible representation matrices of $E$).

We denote by $\mathcal{L}_{A(B)}$ the geometric object consisting of the lattice sites in subsystem $A$ (or in its complement $B)$. Consider the subgroup $\tilde{\mathcal{G}}$ of lattice symmetry operations that can be written as

[TABLE]

where $R^{(A)}$ ($R^{(B)}$) acts only on sites in $A$ ($B$). By

[TABLE]

we denote the corresponding operator that acts onto wave functions in second quantization rather than lattice sites (cf. Eq. (29) below). Then the general invariance of the global $\rho$ under all elements $R$ of the point group $\mathcal{G}$

[TABLE]

implies for the subgroup $\tilde{\mathcal{G}}$ of elements $R^{\prime}$ which can be written in the specific form Eq. (18) that

[TABLE]

In the last step the basis independence of the trace operation and the definition of the reduced density matrix was used. Thus, $\rho_{A}$ is invariant under all joint lattice symmetries $R^{\prime}$ of $\mathcal{L}_{A}$ and $\mathcal{L}_{B}$ that map each subset of lattice sites separately back onto itself according to Eq. (18). If either $\mathcal{L}_{A}$ or $\mathcal{L}_{B}$ has reduced symmetry (e.g. a square plaquette embedded in a rectangular system or a rectangular plaquette inside a square system), then only the largest common symmetry subgroup $\tilde{\mathcal{G}}$ of the point group $\mathcal{G}$ is inherited by $\rho_{A}$. An illustrative example is shown in the inset of Fig. 3 where the plaquette $\mathcal{L}_{A}$ possesses the full symmetry of the square, but due to its location at the corner of a system with open boundary conditions the complement lattice $\mathcal{L}_{B}$ is only invariant under $C_{2x\bar{y}}$, which is reflected in the symmetries of the diagonal elements of $\rho_{A}$ (see main panel of Fig. 3).

Particle number sectors with $N_{\uparrow}=N_{\downarrow}$ can be further decomposed according whether the states are even or odd under spin the inversion symmetry $S$ (see The symmetry group of these sectors is $D_{4h}=D_{4}\times S$, allowing a finer symmetry labelling (see Appendix C).

Once all basis vectors $|\phi_{i\lambda}^{(n)}\rangle$ have been constructed, the reduced density matrix can be transformed from the occupation to the representation basis via

[TABLE]

with the transformation matrix $S_{\alpha,(n,i,\lambda)}=\langle\alpha|\phi_{i\lambda}^{(n)}\rangle$. Here, $|\phi_{i\lambda}^{(n)}\rangle$ is the basis state corresponding to the $\lambda$-th copy of $i$-th row of the $n$-th irreducible representation. The group theoretical techniques for constructing these basis vectors are detailed in Appendix B.

As can be seen from Fig. 2, when written in the irreducible representation basis, the reduced density matrix $\rho_{A}$ does not have non-vanishing matrix elements between states of different symmetry and acquires a block diagonal form. More importantly, we are in a position to attach symmetry labels to the eigenvalues of $\rho_{A}$.

IV Error bars

Due to statistical fluctuations $\rho_{A}$ cannot be perfectly Hermitian, however, the deviations from Hermiticity, $\Delta_{H}=\rho_{A}-\rho_{A}^{\dagger}$, are smaller than the error bars of the corresponding off-diagonal elements, and $\rho_{A}$ is found to be normalized, $\text{Tr}(\rho_{A})=1\pm\varepsilon$, with an inaccuracy ranging from $\varepsilon\lesssim 10^{-5}$ ($0\leq U\leq 4$) to $\varepsilon\approx(1-3)\times 10^{-2}$ (large $U$, $4\leq U\leq 10$) for $\beta t=4$. For low temperatures the inaccuracy is slightly larger ($\epsilon\approx 3\times 10^{-2}$) for all values of $U$ (see insets in Fig. 6 below). Furthermore, $\rho_{A}$ is positive semi-definite within statistical uncertainty, as required for a valid density matrix. Error bars have been obtained with the bootstrap method, in which the matrix diagonalization is repeated $\sim 10^{3}$ times, each time adding Gaussian noise with a standard deviation of the size of the deviation from Hermiticity $|\langle\beta|\Delta_{H}|\alpha\rangle|$ to each matrix element $\langle\beta|\rho_{A}|\alpha\rangle$. The well-resolved symmetry-related degeneracies seen in Figs. 6 (a) and (b) below indicate that this error analysis is sound 444Due to the normalization all elements of $\rho_{A}$ have correlations among them so that adding Gaussian noise idependently to all elements cannot be entirely correct..

V QST for a plaquette in the Hubbard model

A simple argument Scalapino and Trugman (1996) for an isolated plaquette already shows how local antiferromagnetic correlations favour $d$-wave pairing correlations, namely the operator that connects the antiferromagnetic $4$-particle ground state $|4\rangle$ to the two-hole ground state $|2\rangle$ must have $d_{x^{2}-y^{2}}$ symmetry. The matrix element

[TABLE]

is large when the pairing operator $\Delta_{d}=(c_{3,\uparrow}c_{2,\downarrow}-c_{3,\uparrow}c_{4,\downarrow}+\cdots)$ has the sign structure for $d$-wave symmetry. On the other hand $\langle 2|\Delta_{s}|4\rangle=0$ for an $s$-wave pairing operator. The following sections investigate the quantum state on a plaquette of the Hubbard model embedded in a bath of $12\times 12$ sites which are treated numerically exactly.

Probabilities of individual occupation number configurations are shown in Sect. V.1. Sects. V.2 and V.3 present the eigenstates of the plaquette reduced density matrix, with local correlations resolved according to symmetry sectors.

V.1 Diagonal elements of $\rho_{A=\Box}$

Fig. 3 shows the probabilities $P(s)$ of all plaquette configurations $s$ on a plaquette which is located at the corner of a system with open boundary conditions. The integer $s\in\{0,1,\ldots,255\}$ encodes the Fock configuration on a plaquette through its binary representation $[b(s)]=[n_{4}^{\uparrow}n_{3}^{\uparrow}n_{2}^{\uparrow}n_{1}^{\uparrow}n_{4}^{\downarrow}n_{3}^{\downarrow}n_{2}^{\downarrow}n_{1}^{\downarrow}]$ where $n_{i}^{\sigma}$ is the occupation number for spin $\sigma$ at one of the four sites (shown in the upper right inset) and the notation [ ] converts integer codes into bit representations.

The most probable states, the states with integer code [105] and [150], are the two Néel states; they are followed by the 12 other spin-only states which together would span the Hilbert space in a Heisenberg-like description. From the upper left panel of Fig. 3 one can see that for $U=7.2$, $T/t=0.35$ and half filling, on a plaquette there are $\sim 24\%$ Néel states, $\sim 40\%$ spin-only states (excluding the two Néel states), the remaining $\sim 36\%$ are states with charge fluctuations.

The arrangement with the plaquette at the corner, shown in the upper right inset in Fig. 3, does not possess the full symmetry of the square, the only symmetry operations which respect Eq. (18) being $\{E,C_{2x\bar{y}}\}$. This is reflected in asymmetries of the probabilities for plaquette configurations with a single hole (see main panel of Fig. 3 with plaquette configurations drawn next to representative data points): In the presence of a boundary the hole prefers to have many neighbours rather than sit at the boundary which would limit the number of possible hopping processes. Therefore, among the configurations shown, the one with the hole located precisely at the corner has the lowest probability. Note that configurations that are related by the symmetry operation $\{E,C_{2x\bar{y}}\}$ do occur with the same probability.

Fig. 4 shows the probabilities of selected plaquette occupation number states in the repulsive Hubbard model at half filling for low temperature $\beta t=16$ or $24$ (a) and high temperature $\beta t=4$ (b) as a function of the interaction strength $U/t$.

Here, periodic boundary conditions are used so that the full symmetry of the square is preserved. Occupation number states that are related by symmetries are grouped into classes of states, which are labelled by the bitcode $[b(s)]$ of the member with the smallest bitcode within the class. A list of all 34 classes of symmetry-related states with the bitcodes of their representatives can be found in Appendix A.

At low temperature [Fig. 4(a)], two datasets for different temperatures, $\beta t=24$ for $U/t\leq 5$ and $\beta t=16$ for $U/t\geq 6$, have been combined, which is indicated by different colours of the background shading. For $U/t\geq 3$, the most probable states are the two Néel states. The second most important class of states [30] comprises those states with three particles of one spin and one particle of the opposite spin. Taking into account spin-inversion symmetry there are 8 such states and their probabilities have been summed up in Fig. 4. This explains why this class of states appears with higher probability in Fig. 4 than the class of states [60], consisting of only 4 states, namely the spin configurations that are FM aligned in one direction and AFM aligned in the other one. The multiplicities of different classes of states, which are essential for interpreting Fig. 4, are also given in Appendix A.

To illustrate that detailed information can be gleaned from configurational probabilities we consider next in Fig. 4 the classes [27] and [54], which both have a multiplicity of 16 states. The class of states [27], comprising plaquette configurations with a neighbouring doublon-hole pair and the singly-occupied sites in a FM configuration, has smaller probability than states [54] representing a neighbouring doublon-hole pair with the singly-occupied sites arranged in an AFM configuration. Note also that the classes of states [19] and [55] have exactly the same probability due to particle-hole symmetry at half filling.

The plaquette probabilities at high temperature $\beta t=4$ [Fig. 4(b)] are qualitatively very similar to those at low temperature, indicating that local correlations of the low-temperature phase are already well developed at $\beta t=4$.

In Fig. 5(a) the doping dependence of the plaquette probabilities is presented for experimentally relevant inverse temperature $\beta t=4$ and for repulsive interaction $U/t=2,4,$ and $7.2$. Error bars in Fig. 5(a) are deduced from the spread of datapoints within one class of states which should have the same probability due to symmetry.

The selected plaquette configurations have at least one doubly occupied site so that for large Hubbard repulsion they represent the intermediate virtual states through which pairwise and ring exchange interactions in an effective spin Hamiltonian are mediated Takahashi (1977); MacDonald et al. (1988); Delannoy et al. (2005). A typical pathway of hopping processes leading in fourth order perturbation theory to ring exchange interactions is illustrated in Fig. 5(b).

The overall trend is that the amplitude of states with charge fluctuations is reduced with increasing Hubbard repulsion. The amplitude of states with neighbouring doublon-hole pairs decreases as $t/U$, in accordance with second order perturbation theory. Plaquette configurations with a doublon-hole pair on diagonally opposite corners ([53]), occurring as intermediate states in ring-exchange pathways [Fig. 5(b)], have a very small probability at all fillings, which is for large $U/t$ approximately an order of magnitude smaller than that of states with neighbouring doublon-hole pairs ([27],[54]), as to be expected on the basis of fourth order perturbation theory.

Based on Fig. 5(a), very detailed observations regarding correlation effects can be made. For example, a signal for local correlation lies in the differences of probability between similar configurations such as [23] and [53] or [19] and [25], which would be equally likely, if the singly and doubly occupied sites were placed on the lattice randomly, with say probability $p_{d}$ for a doubly occupied site and $p_{s}$ for a singly occupied one. It is important to note that [19] has multiplicity 16, whereas [25] has multiplicity 8 (see Appendix A). Thus, the observation that the classes [19] and [25] have almost the same probabilities in Fig. 5(a) indicates that the probability per individual configuration of states [25] is kinetically enhanced compared to states from [19], since the former allow for more hopping processes on the plaquette. Note also in Fig. 5(a) at $U/t=7.2$ the discontinuous jump at quarter filling $\langle n\rangle\approx 0.5$.

V.2 Off-diagonal elements of $\rho_{A=\Box}$

With the knowledge of all off-diagonal elements of the reduced density matrix, we can compute the entanglement spectrum and resolve it according to symmetry sectors.

Fig. 6 shows the eigenvalue spectrum of the plaquette reduced density matrix $\rho_{A=\Box}$ at half filling as a function of Hubbard interaction for high (a) and low (b) temperature. The 16 most important eigenstates of $\rho_{A=\Box}$ are labelled by their symmetry sectors $(N_{A,\uparrow},N_{A,\downarrow})M$ where $(N_{A,\uparrow},N_{A,\downarrow})$ is the particle number sector and $M$ is the Mulliken symbol describing the symmetry of the state under the operations of the symmetry group $D_{4h}$ which combines the symmetries of the square and spin inversion symmetry. The data for high temperature in Fig. 6(a) is presented for the repulsive Hubbard model, while for the low temperature data in Fig. 6(b) we have chosen the language of the attractive Hubbard model. At half filling, both models are exactly equivalent; the eigenstates of $\rho_{A=\Box}$ for the repulsive Hubbard model are related to those of the attractive model by the spin-down particle-hole transformation

[TABLE]

where $(-1)^{{\bf i}}\equiv(-1)^{i_{x}+i_{y}}$ is a staggered phase factor for one of the two sublattices of the square lattice. We point out that equivalent states in both models behave differently under symmetry operations of the point group, which is reflected in their Mulliken symbols. There is a one-to-one correspondence between the Mulliken symbols of equivalent states in the repulsive and attractive model, for example states with symmetry $(2,2)A_{2,u}$ in the language of the positive-$U$ model correspond to states with symmetry $(2,2)B_{2,g}$ in the negative-$U$ model, etc. We stress that in the presentation of Fig. 6(a) and (b) same colours do not necessarily imply that the states are equivalent in the two models 555 The colours are based on the sequence of energy levels at $|U|/t=4$..

The insets in Fig. 6(a) and (b) show the normalization $\text{Tr}\left(\rho_{A=\Box}\right)$ and the purity $\text{Tr}\left(\rho_{A=\Box}^{2}\right)$ of the plaquette reduced density matrix. The purity increases with interaction strength, signalling that the state becomes closer to the $T=0$ limiting case of a product state of non-entangled plaquettes. For non-interacting fermions ($U=0$) the plaquette is most mixed with its environment; then there is a pronounced upturn in $\text{Tr}\left(\rho_{A=\Box}^{2}\right)$ around $U=5$, which levels off for $U\rightarrow 10$. The purity at $\beta t=4$ (Fig. 6(a)) is by a factor of two smaller than at $\beta t=24$ (Fig. 6(b)), as to be expected due to thermal entropy. $\rho_{A=\Box}$ is normalized within error bars, which, however, increase for larger $U$ and lower temperature.

Focusing first on results for high temperature Fig. 6(a), we observe that at $U=0$ the 16 most important eigenstates are all degenerate and clearly separated from the remaining eigenstates with lower weight. As the interactions are switched on, the multiplet splits into two singlet-triplet type sequences, namely the four states $|s^{(1)}\rangle\in(2,2)B_{1,g}$ and $\left\{|t^{(1)}_{0}\rangle\in(2,2)A_{2,u},\,|t^{(1)}_{-}\rangle\in(1,3)A_{2},\,|t^{(1)}_{+}\rangle\in(3,1)A_{2}\right\}$, and the four states $|s^{(2)}\rangle\in(2,2)A_{1,g}$ and $\left\{|t^{(2)}_{0}\rangle\in(2,2)B_{2,g},\,|t^{(2)}_{-}\rangle\in(1,1)A_{1,g},\,|t^{(2)}_{+}\rangle\in(3,3)A_{1,g}\right\}$, and into a degenerate octet of eigenstates with $p$-wave symmetry (labelled by the irreducible representation label $E$).

The states with $p$-wave symmetry are both spin and pseudospin doublets, which together with the fact that their irreducible representation $E$ is two-dimensional explains their eightfold degeneracy Schumann (2002). The degeneracy will be lifted either by applying an external magnetic field or by shifting the chemical potential away from the half filling point $\mu=\frac{U}{2}$. It is remarkable how well the degeneracy of the octet is preserved in the Monte Carlo data of Fig. 6(a) and (b).

In comparison with the grand canonical eigensystem of an isolated 4-site Hubbard model, which is worked out analytically in Ref. Schumann (2002), the succession of energy levels appears changed in Fig. 6(a) and (b) in that the “entanglement energy” of the octet is lower (i.e. it has a heigher weight in the thermal state) than that of the low-weight singlet-triplet-structure $\{|s^{(2)}\rangle,|t^{(2)}_{\pm,0}\rangle\}$ for $|U|<6$. This shows that the method presented here can resolve fine differences between the energy spectrum of an isolated plaquette Schumann (2002) and of a plaquette embedded in a much larger system. Note that it is not simply the difference in temperature between the ground state spectrum of Ref. Schumann (2002) and our “entanglement energy“ spectrum, which could explain this discrepancy, since it occurs also at $\beta=24$ in Fig. 6(b) and since temperature cannot change the relative order of the statistical weights.

For $U>4$, the degeneracy of the high-weight triplet $|t^{(1)}_{\pm,0}\rangle$ is lifted, and concomitantly the purity of the reduced density matrix increases (see inset Fig. 6 (a)). The fact that the octet of eigenstates remains perfectly degenerate in the interval $4\leq U\leq 6$ supports the picture that the lifting of the degeneracy in $|t^{(1)}_{\pm,0}\rangle$ is not merely an artifact of larger error bars. The low-weight singlet state $|s^{(2)}\rangle\in(2,2)A_{1,g}$ shows non-monotonic behaviour as a function of $U$ and its weight appears to increase again for $U>6$. However, the error bars are too large to draw any conclusions.

In Fig. 6(b) entanglement spectra for two low temperatures, $\beta t=24$ (for $|U|/t\leq 5$) and $\beta t=16$ (for $|U|/t\geq 6$), are combined. Due to issues of ergodicity at large Hubbard interactions it was not possible to reach lower temperatures for $|U|/t\geq 6$. The eigenvalue spectra $\{\lambda_{n}\}$ have a qualitatively similar dependence on $U$ both for high and low temperature. In both cases, multiplets of degenerate states that exist in the sector with charge fluctuations for small $U$ mix for large $|U|/t$ and give rise to a broad structureless ”band“ of small eigenvalues, from which a low-lying singlet state seems to separate off. However, large error bars prevent a conclusive statement.

When discussing now to the low-temperature spectrum of $\rho_{A=\Box}$ displayed in Fig. 6(b), we refer to the Mulliken symbols shown next to Fig. 6(b). It must be stessed again that equivalent states in the repulsive and attractive Hubbard model, i.e. states related by the spin-down particle-hole transformation Eq. (24), are labelled by different Mulliken term symbols. For the positive-$U$ Hubbard model, the singlet-triplet type structures are comprised of the four states $|s^{(1)}\rangle\in(2,2)A_{1,g}$ and $\{|t_{0}^{(1)}\rangle\in(2,2)B_{2,g},\,|t_{+}^{(1)}\rangle\in(3,3)A_{1,g},\,|t_{-}^{(1)}\rangle\in(1,1)A_{1,g}\}$ and the four states $|s^{(2)}\rangle\in(2,2)B_{1,g}$ and $\{|t_{0}^{(2)}\rangle\in(2,2)A_{2,u},\,|t_{+}^{(2)}\rangle\in(3,1)A_{2},\,|t_{-}^{(2)}\rangle\in(1,3)A_{2}\}$. In between the two singlet-triplet structures there is again an octet of degenerate states. A remarkable difference to the high-temperature spectrum is that the degeneracy of the high-weight triplet is lifted at smaller $|U|/t$, namely at $|U|/t=2$ for $\beta t=24$ compared to $U/t=4$ for $\beta t=4$. The weight of $|t_{0}^{(1)}\rangle$ increases rapidly, the subsequent decrease at $|U|/t=6$ must be attributed to the change in temperature from $\beta t=24$ to $\beta t=16$ when changing from the dataset with $\beta t=24$ for $|U|/t\leq 5$ to the dataset with $\beta t=16$ for $|U|/t\geq 6$.

V.3 Doping dependence

Next, we turn to the doping dependence of the plaquette entanglement spectrum, which is displayed in Fig. 7 for the repulsive Hubbard model at $U/t=6$ and high temperature $\beta t=4$. The DQMC algorithm suffers from a sign problem Loh Jr. and Gubernatis (1992) in the repulsive Hubbard model when particle-hole symmetry is broken by tuning the chemical potential away from the half-filling point $\mu=\frac{U}{2}$. Nevertheless, at relatively high temperature around $\beta t=4$ simulations are still possible due to a mild sign problem Iglovikov et al. (2015) that can be offset by acquiring more statistics in longer Monte Carlo runs. Fig. 7 shows the eigenvalue spectrum $\{\lambda_{n}\}_{n=1}^{4^{N_{s}}}$ of the plaquette $(N_{s}=4)$ reduced density matrix $\rho_{A=\Box}$ with coloured stripes indicating blocks of fixed particle number $(N_{\uparrow},N_{\downarrow})$.

In the following, the most important eigenstates from the respective symmetry multiplets are listed explicitly. At half filling, the leading eigenstate of the plaquette reduced density matrix has $d_{x^{2}-y^{2}}$-wave ($B_{1g}$) symmetry, whereas around quarter filling, $\langle n\rangle\approx 0.5$, it has $s$-wave ($A_{1g}$) symmetry (see first and fourth row of Fig. 7). The invariant subspace labelled by $(2,2)B_{1g}$ consists of four states and $(2,2)A_{2u}$ consists of three states (see Tab. 3), but it turns out that a single state from each symmetry multiplet has by far the largest coefficient, namely:

[TABLE]

Dots in Eq. (25a) indicate states from the same symmetry multiplet without double occupancy but with much smaller weight. Thus, the leading eigenvectors of $\rho_{A=\Box}$ in the particle number sector $(2,2)$ at half filling are the symmetric and antisymmetric combinations of the two Néel states, as to be expected. Note that for $U\gtrsim 6$ an additional $p$-wave doublet $(E_{u}^{x},E_{u}^{y})$ appears in the particle number sector $(2,2)$ (see first row of Fig. 7). The leading eigenvectors in the particle number sector $(1,1)$ at quarter filling are

[TABLE]

i.e. two holes in diagonally opposite corners, which maximizes their kinetic energy, and a set of states with smaller weight with two neighbouring holes, which have $p$-wave symmetry (labelled by $E_{u(g)}^{x(y)}$ in the fourth row in Fig. 7). The character of the leading eigenstates in the particle number sectors $(2,2)$ and $(1,1)$ as described by Eqs. (25a) and (26) hardly changes with doping. The leading eigenstate from the symmetry sector $(1,3)A_{2}$

[TABLE]

and its spin-reversed counterpart from the symmetry sector $(3,1)A_{2}$ are degenerate with the leading eigenstate from $(2,2)A_{2,u}$ (Eq. (25b)) up to $U/t\lesssim 4$ (at $\beta t=4$). Finally, we list the state

[TABLE]

with spins that are ferromagnetically aligned along one coordinate axis and antiferromagnetically in the other direction and which separates off from the “band” of low-lying states for large $|U|$ (see Fig. 6(a)).

VI Conclusion and outlook

We have provided proof of principle calculations that it is possible in an equilibrium DQMC simulation to obtain the full reduced density matrix of a small subsystem embedded in a much larger system that can be interpreted as the exact correlated bath. The sequence of “entanglement energy” levels of the embedded plaquette is shown to differ from the sequence of levels for an isolated plaquette Schumann (2002).

The calculated configurational probabilities allow detailed benchmarking of current fermionic quantum gas microscope experiments.

The possiblity of computing the full quantum state of a subsystem is a unique feature of the DQMC framework, which is based on the free fermion decomposition Grover (2013). Due to the factorization of the Monte Carlo weight into a spin-$\uparrow$ and spin-$\downarrow$ part, the computational cost for obtaining all elements of $\rho_{A}$ scales like $2\times(2^{N_{s}}\times 2^{N_{s}})$, rather than $4^{N_{s}}\times 4^{N_{s}}$, in one Monte Carlo sample. However, the storage requirement for all elements is $4^{N_{s}}\times 4^{N_{s}}$, which is forbidding for e.g. $N_{s}=9$. If individual (particle number and point group) symmetry sectors are targeted by performing the transformation Eq. (22) in every Monte Carlo step rather than computing all elements of $\rho_{A}$ in the occupation number basis, the number of non-vanishing matrix elements of $\rho_{A}$ that need to be kept for Monte Carlo averaging can be reduced to $a_{n}\times a_{n}$, where $a_{n}$ is the number of copies of the $n$-th irreducible representation that appear in the decomposition of a given particle number sector (see Eq. (31)). The additional cost of the basis transformation Eq. (22) in every Monte Carlo step can be compensated by a finer granularity of the parallelization. This would give access to $3\times 3$ subsystems that can already capture the effect of next nearest neighbour hopping $t^{\prime}$, which needs to be included to describe qualitatively the electronic band structure of cuprates. At least for temperatures, where DQMC simulations are still possible in spite of the sign problem Iglovikov et al. (2015), one may thus hope to gain some insight into the role of local correlations in the high-temperature phase of a prototypical model for high-$T_{c}$ superconductors, whose pseudogap regime and anomalous normal state commonly referred to as “strange” metal phase is still poorly understood Keimer et al. (2015).

A model for which quantum state tomography on a plaquette is particularly meaningful is the plaquette Hubbard model studied in Refs. Rey et al. (2009) and Ying et al. (2014), which interpolates between isolated plaquettes and a uniform square lattice taking the interplaquette hopping as a tunable parameter. The presented approach may also prove useful for computing dynamical properties such as the spectral function or optical conductivity without the need for analytical continuation of imaginary-time correlation functions, in some form of cluster approximation (see e.g. Ref. Sénéchal et al. (2000)), albeit with an exact correlated bath.

Acknowledgments

The author thanks Lei Wang for helpful discussions and for pointing out Refs. Udagawa and Motome (2010, 2015). The numerical simulations were performed on JURECA, Jülich Supercomputing Center. Support by the International Young Scientist Fellowship of Institute of Physics, Chinese Academy of Sciences under the Grant No. 2018004 is acknowledged.

Appendix A List of symmetry-related classes of states for a single plaquette

Each of the 256 plaquette states is labelled by an integer $x$ between 0 and 255 with the convention that its binary representation $[x]$ corresponds to the occupation numbers on the plaquette; the four least significant bits denote occupation numbers for spin-$\downarrow$ (see Fig. 8).

The leftmost column of Fig. 8 shows microscpic configurations of spin-$\uparrow$ and spin-$\downarrow$ particles on the plaquette. States that are related by a symmetry operation of the point group $D_{4}$ of the square or by spin-inversion are grouped into classes of symmetry-related states. The 34 classes are listed with a representative spin configuration for each class, its bitcoded label and the number of symmetry-related states in the class (“multiplicity”). We choose as a representative from each class the state with the smallest bitcoded label. Special sets of states are the two Néel states (labelled by $[105]$), the 16 states in the spin-only subspace (with representatives $[15]$,$[30]$,$[60]$ and $[105]$), the states with neighbouring doublon-hole pairs in a spin-only background ($[27]$, $[54]$) and the states with diagonally-opposite doublon-hole pairs in a spin-only background ($[23]$,$[53]$).

Appendix B Group theoretic techniques: projection operator method

The relation between transformation operators acting in the many-body Hilbert space and the symmetry operations acting on coordinates is provided by Wigner’s convention (Tinkham, 1964, Chapt. 3)

[TABLE]

where $\hat{P}_{R}$ is the operator acting onto wave functions in second quantization while the symmetry operator $R$ acts onto indices of creation and annihilation operators.

The fermion ordering in the definition of the states is chosen such that site indices of creation operators increase from right to left and creation operators for $\uparrow$-particles are to the left of operators for $\downarrow$-particles, e.g. $|\uparrow_{4}\,\downarrow_{3}\,\uparrow_{2}\,\downarrow_{1}\,\rangle\equiv c_{4\uparrow}^{\dagger}c_{2\uparrow}^{\dagger}c_{3\downarrow}^{\dagger}c_{1\downarrow}^{\dagger}|\text{vac}\rangle$. As an illustration of Eq. (29) and of the action of the symmetry operators $\hat{P}_{R}$ on the many-body Hilbert space, consider the example

[TABLE]

which shows how the matrix elements of $\hat{P}_{R}$ can be constructed. In the particle number block $(N_{\uparrow},N_{\downarrow})$, the operators $\hat{P}_{R}$ are permutation matrices of size $\text{dim}(N_{\uparrow},N_{\downarrow})\times\text{dim}(N_{\uparrow},N_{\downarrow})$, where $\text{dim}(N_{\uparrow},N_{\downarrow})=\begin{pmatrix}N_{s}\\ N_{\uparrow}\end{pmatrix}\cdot\begin{pmatrix}N_{s}\\ N_{\downarrow}\end{pmatrix}$, with an additional sign structure coming from the fermionic exchanges. Since the spatial symmetry operations do not affect the spin states, it is convenient to write $\hat{P}_{R}$ as the tensor product $\hat{P}_{R}=\hat{P}_{R,\uparrow}\otimes\hat{P}_{R,\downarrow}$ with $\hat{P}_{R,\sigma}$ acting only on creation operators of spin $\sigma$. Having obtained a matrix representation of the symmetry operators $\hat{P}_{R}$ on the particle number sector $(N_{\uparrow},N_{\downarrow})$, we can decompose this subspace of Hilbert space further into the irreducible invariant subspaces of $D_{4}$ via the projection operator technique Tinkham (1964) (see also Kuns et al. (2011) for a detailed discussion). In the decomposition of a reducible representation the $n$-th irreducible representation occurs $a_{n}$ times, given by Tinkham (1964)

[TABLE]

where $\chi^{(n)}(R)$ is the character of the group element $R$ in the $n$-th irreducible representation and $\chi(R)\equiv\text{Tr}\left(P_{R}\right)=\sum_{i}\left[P_{R}\right]_{ii}$ is the character of $R$ in the reducible matrix representation. Applying the formula (31) to each particle number sector $(N_{\uparrow},N_{\downarrow})$ of a square plaquette we obtain the group structure presented in Tab. 2.

We wish to decompose the particle number sector $(N_{\uparrow},N_{\downarrow})$ into blocks of states such that the application of a lattice symmetry operation to a state mixes only states within the same block.

Let $|\phi_{i\lambda}^{(n)}\rangle$ denote a normalized basis state that transforms according to the $\lambda$-th copy of the $i$-th row in the $n$-th irreducible representation. Then each occupation number state $|\alpha\rangle$ can be expanded as

[TABLE]

where $c$ is the number of irreducible representations, which is equal to the number of conjugacy classes Tinkham (1964) (here, for $D_{4}$, c=5), $l_{n}$ is the dimension of the $n$-th irreducible representation, and $\lambda$ labels the $a_{n}$ different copies of the $n$-th irreducible representation.

The symmetry transfer operator is defined as Tinkham (1964)

[TABLE]

where $\Gamma^{(n)}(R)$ is the matrix representation of the group element $R$ in the $n$-th irreducible representation and the sum runs over all group elements. $\mathcal{P}_{ii}^{(n)}$ acts as a projector onto the $i$-th row of the $n$-th irreducible representation, while $\mathcal{P}_{ij}^{(n)}$ transfers the $i$-th row into the $j$-th row according to Tinkham (1964)

[TABLE]

and

[TABLE]

Note that if there are several copies $\lambda$ of the same irreducible representation $n$, then the projection operator $\mathcal{P}_{ii}^{(n)}$ applied to a basis state $|\alpha\rangle$ will return a basis state for only a single copy $\lambda(\alpha)$:

[TABLE]

By letting $\mathcal{P}_{ii}^{(n)}$ act onto each state $|\alpha\rangle$ of the particle number sector $(N_{\uparrow},N_{\downarrow})$ and collecting all non-zero states that are linearly independent, all copies of the $n$-th irreducible representation are generated. In this projection method, it may happen that the same basis state (up to a global phase) is generated multiple times. Thus, the basis vectors of all irreducible representations can be constructed and combined into a unitary matrix $S_{\alpha,(n,i,\lambda)}=\langle\alpha|\phi_{i\lambda}^{(n)}\rangle$ which transforms the reduced density matrix from the occupation number basis to the representation basis according to Eq. 22 of the main text.

Appendix C Spin inversion symmetry

For a finer symmetry labelling it is useful to implement the spin inversion symmetry $\mathcal{S}=\{\mathcal{E},\sigma_{h}\}$ where $\sigma_{h}=\prod_{i\in A}\sigma_{i}^{x}$ flips all spins on subsystem $A$. All symmetry operations of the lattice symmetry group $D_{4}$ commute with the spin inversion operation since they act onto different degrees of freedom (site indices of creation operators on the one hand and spin indices on the other hand). Therefore, we can form the direct-product group $D_{4h}=D_{4}\times\mathcal{S}$ with 16 group elements, the original 8 from $D_{4}$, each multiplied by the identity or by spin inversion, and organize the states into the irreducible invariant subspaces of $D_{4h}$. In order to apply the projection operator method for generating the irreducible basis states, one needs to know the irreducible representation matrices of $D_{4h}$ (see Eq.(33)). It can be shown Tinkham (1964) that the direct product of two irreducible representations forms an irreducible representation of the direct product group.

If the Hilbert space is first decomposed into subspaces of fixed particle number $H=\prod_{\oplus\\ N_{\uparrow},N_{\downarrow}=0}^{N_{s}}H_{(N_{\uparrow},N_{\downarrow})}$, then spin inversion symmetry $\mathcal{S}$ can only be used for further block diagonalization inside subspaces with equal spin populations since it is obviously not possible to construct eigenstates of $\sigma_{h}$ that lie only in $H_{(N_{\uparrow},N_{\downarrow})}$ whenever $N_{\uparrow}\neq N_{\downarrow}$. Including spin inversion leads to the finer group structure of the subspaces with $N_{\uparrow}=N_{\downarrow}=1$ and $N_{\uparrow}=N_{\downarrow}=2$ shown in Tab. 3 where the additional label $g$ ($u$) indicates whether the basis function is even (odd) under spin inversion.

Appendix D Pseudocode for computing $\langle\beta|\rho_{A}|\alpha\rangle_{\{{\bf s}\}}$

The following code listing describes how to compute a matrix element of the reduced density matrix $\rho_{A}$ for a free fermion system in the external potential of auxiliary fields $\langle{\bf s}\rangle$ given the single particle Green’s function on subsystem $A$ as input. The main task consists in collecting the appropriate row and column indices for the submatrices, which appear in the determinant formula of Eq. (9). The symbols $\mathcal{I}_{+}$ and $\mathcal{I}_{-}$ and the meaning of the abbreviations ”occ“ and ”ua“ are defined in the main text. The phase factors in line 23 are those from Eqs. (14) and (15) in the main text. Note that the computation of a many-body reduced density matrix from a single-particle Green’s function obviously constitutes a blow-up of redundant information.