Full Counting Statistics for Interacting Fermions with Determinantal Quantum Monte Carlo Simulations

TL;DR

This paper introduces a method to compute the full probability distribution of quadratic observables in the Fermi-Hubbard model using determinantal quantum Monte Carlo, aiding the analysis of cold atom experiments and magnetic properties.

Contribution

The paper develops a novel approach for calculating full counting statistics in interacting fermion systems within the determinantal quantum Monte Carlo framework.

Findings

Full counting statistics reveal information on preformed pairs.

Excellent agreement with experimental measurements of staggered magnetization.

Current experiments can distinguish between Hubbard and Heisenberg models.

Abstract

We present a method for computing the full probability distribution function of quadratic observables for the Fermi-Hubbard model within the framework of determinantal quantum Monte Carlo. Especially, in cold atoms experiments with single site resolution, such full counting statistics can be obtained from repeated projective measurements. We demonstrate, that the full counting statistics can provide important information on the size of preformed pairs. Furthermore, we compute the full counting statistics of the staggered magnetization in the repulsive Hubbard model at half filling and find excellent agreement with recent experimental results. We show that current experiments are capable of probing the difference between the Hubbard model and the limiting Heisenberg model.

Equations

$$\chi ({\varphi }^{\uparrow },{\varphi }^{\downarrow })=⟨e^{i{\varphi }^{\uparrow }{\hat{N}}_A^{\uparrow }+i{\varphi }^{\downarrow }{\hat{N}}_A^{\downarrow }}⟩=\text{Tr}\left({\rho }_A{e}^{i{\sum }_{\sigma }{\varphi }^{\sigma }{\hat{N}}_A^{\sigma }}\right).$$

$$P(N_A^{\uparrow },N_A^{\downarrow })=\sum _{n,m=0}^{N_s}\frac{e^{-i{\varphi }_n^{\uparrow }{N}_A^{\uparrow }-i{\varphi }_m^{\downarrow }{N}_A^{\downarrow }}}{(N_s+1{)}^2}\chi ({\varphi }_n^{\uparrow },{\varphi }_m^{\downarrow })$$

$$P_N(N_A)=\sum _{N_A^{\downarrow }=0}^{N_s}P(N_A-N_A^{\downarrow },N_A^{\downarrow })$$

$$⟨\hat{O}⟩=\text{Tr}\left(\hat{O}{e}^{-\beta \hat{H}}\right)/\text{Tr}\left(e^{-\beta \hat{H}}\right)=\frac{1}{Z}\sum _{\{\mathbf{s}\}}{w}_{\{\mathbf{s}\}}{O}_{\{\mathbf{s}\}},$$

$$w_{\{\mathbf{s}\}}=\det {\left({\mathcal{G}}_{\{\mathbf{s}\}}^{\uparrow }\right)}^{-1}{\left({\mathcal{G}}_{\{\mathbf{s}\}}^{\downarrow }\right)}^{-1},$$

$$\chi ({\varphi }^{\uparrow },{\varphi }^{\downarrow })=\sum _{\{\mathbf{s}\}}{w}_{\{\mathbf{s}\}}{\chi }_{\{\mathbf{s}\}}^{\uparrow }\left({\varphi }^{\uparrow }\right){\chi }_{\{\mathbf{s}\}}^{\downarrow }\left({\varphi }^{\downarrow }\right)$$

$${\rho }_{\{\mathbf{s}\},A}^{\sigma }={\mathcal{K}}_{\{\mathbf{s}\},A}^{\sigma }{e}^{-H_A^{\sigma }}$$

$$H_A^{\sigma }=-\sum _{i,j\in A}{\hat{c}}_{i,\sigma }^{†}\log {\left({\left[G_{\{\mathbf{s}\},A}^{\sigma }\right]}^{-1}-1\right)}^{ij}{\hat{c}}_{j,\sigma }$$

$${\chi }_{\{\mathbf{s}\}}^{\sigma }(\varphi )=\prod _{\alpha =1}^{N_s}(1-{\lambda }_{\alpha }^{\sigma })\prod _{\alpha =1}^{N_s}\left(1+e^{-{ϵ}_{\alpha }^{\sigma }+i\varphi }\right),$$

$${\chi }_{\{\mathbf{s}\}}^{\sigma }(\varphi )=\prod _{\alpha =1}^{N_s}(1+(e^{i\varphi }-1){\lambda }_{\alpha }^{\sigma }),$$

$$\mathcal{H}=-t\sum _{⟨\mathbf{r},{\mathbf{r}}^{\prime }⟩,\sigma }\left({\hat{c}}_{\mathbf{r},\sigma }^{†}{\hat{c}}_{{\mathbf{r}}^{\prime },\sigma }+h.c.\right)+U\sum _{\mathbf{r}}{\hat{n}}_{\mathbf{r},\uparrow }{\hat{n}}_{\mathbf{r},\downarrow }.$$

$$P_{\text{odd}}(L_A)=\frac{1}{2}-\frac{1}{2}(1-2p_{\text{deloc}}{)}^{4L_A\cdot ⟨n⟩}$$

$$\begin{gathered} {\chi }_{\{\mathbf{s}\}}^{st}(\varphi )={\text{Tr}}_{\uparrow }\left(e^{-H_A^{\uparrow }}{e}^{i\varphi {\sum }_{\text{i}\in A}(-1{)}^{\mathbf{i}}{\hat{n}}_{i,\uparrow }}\right) \\ {\text{Tr}}_{\downarrow }\left(e^{-H_A^{\downarrow }}{e}^{i\varphi {\sum }_{i\in A}(-1{)}^{\mathbf{i}}{\hat{n}}_{i,\downarrow }}\right). \end{gathered}$$

$${\left[{\tilde{H}}_A^{\sigma }\right]}_{i,j}={\left[H_A^{\sigma }\right]}_{i,j}+i\varphi (-1{)}^{\mathbf{i}}{\delta }_{i,j}$$

$$H_A^{\sigma }=-\sum _{i,j\in A}{\hat{c}}_{i,\sigma }^{†}\log {\left({\left[G_{\{\mathbf{s}\},A}^{\sigma }\right]}^{-1}-1\right)}^{ij}{\hat{c}}_{j,\sigma },$$

$${\chi }_{\{\mathbf{s}\}}^{st}(\varphi )=\prod _{\sigma =\uparrow ,\downarrow }\left[\prod _{\alpha =1}^{N_s}(1-{\lambda }_{\alpha }^{\sigma })\left(1+e^{-{\widetilde{ϵ}}_{\alpha }^{\sigma }(\varphi )}\right)\right].$$

$${\left[H_A^{\sigma }\right]}_{i,j}=-{\left[\log \left(G_{\{\mathbf{s}\},A}^{\sigma }(1-G_{\{\mathbf{s}\},A}^{\sigma }{)}^{-1}\right)\right]}_{i,j}.$$

$$\begin{gathered} {\mathcal{H}}_{U<0}=\sum _{\mathbf{k},\sigma }({\epsilon }_{\mathbf{k}}-\mu ){\hat{c}}_{\mathbf{k},\sigma }^{†}{\hat{c}}_{\mathbf{k},\sigma } \\ -\frac{|U|}{N_{\text{sites}}}\sum _{\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{q}}{\hat{c}}_{\mathbf{k}+\mathbf{q},\uparrow }^{†}{\hat{c}}_{-\mathbf{k}+\mathbf{q},\downarrow }^{†}{\hat{c}}_{{\mathbf{k}}^{\prime }+\mathbf{q},\uparrow }{\hat{c}}_{-{\mathbf{k}}^{\prime }+\mathbf{q},\downarrow }, \end{gathered}$$

$$v_{\mathbf{k}}^2=1-u_{\mathbf{k}}^2=\frac{1}{2}\left(1-\frac{{\epsilon }_{\mathbf{k}}^{(\text{HF})}-\mu (n,U)}{E_{\mathbf{k}}}\right),$$

$$|\text{BCS}⟩=\prod _{\mathbf{k}}\left(u_{\mathbf{k}}(\mu ,\Delta )+v_{\mathbf{k}}(\mu ,\Delta ){\hat{c}}_{\mathbf{k},\uparrow }^{†}{\hat{c}}_{-\mathbf{k},\downarrow }^{†}\right)|\text{vac}⟩,$$

$$\frac{1}{|U|}=\frac{1}{N_{\text{sites}}}\sum _{\mathbf{k}\in \text{1st BZ}}\frac{1}{2\sqrt{({\epsilon }_{\vec{k}}-\mu -|U|\frac{n}{2}{)}^2+{\Delta }^2}},$$

$$n=1-\frac{1}{N_{\text{sites}}}\sum _{\mathbf{k}\in \text{1st BZ}}\frac{{\epsilon }_{\vec{k}}-\mu -|U|\frac{n}{2}}{\sqrt{({\epsilon }_{\vec{k}}-\mu -|U|\frac{n}{2}{)}^2+{\Delta }^2}}.$$

$${\chi }^{(N)}(\varphi )=⟨e^{i\varphi {\hat{N}}_A}⟩=⟨\prod _{i\in A}\prod _{\sigma =\uparrow ,\downarrow }{e}^{i\varphi (1-{\hat{n}}_{i,\sigma })}⟩=⟨\prod _{i\in A}\prod _{\sigma =\uparrow ,\downarrow }{E}_{i\sigma }{F}_{i\sigma }(\varphi )⟩,$$

$${\chi }^{(N)}(\varphi )$$

$$=\prod _{k=1}^{2N_s}({\mu }_k+(1-{\mu }_k)e^{i\varphi }).$$

$$M=\left(\begin{array}{cc}⟨{\hat{c}}_{i\uparrow }^{†}{\hat{c}}_{j\uparrow }{⟩}_{i,j\in A}& ⟨{\hat{c}}_{i\uparrow }{\hat{c}}_{j\downarrow }{⟩}_{i,j\in A}\\ ⟨{\hat{c}}_{i\downarrow }{\hat{c}}_{j\uparrow }{⟩}_{i,j\in A}& ⟨{\hat{c}}_{i\downarrow }^{†}{\hat{c}}_{j\downarrow }{⟩}_{i,j\in A}\end{array}\right).$$

$$p_1\equiv p(m_i^{st}=+1)=p(\uparrow ,+)+p(\downarrow ,-)$$

$$p_2\equiv p(m_i^{st}=-1)=p(\uparrow ,-)+p(\downarrow ,+)$$

$$p_3\equiv p(m_i^{st}=0)=p(d)+p(h)=1-p_1-p_2,$$

$$P(M_A^{st}=k-l)=\sum _{k=0}^{N_s}\sum _{l=0}^{N_s-k}{\delta }_{M_A^{st},k-l}{P}_3(X=k,Y=l)$$

$$P_3(X=k,Y=l)=\frac{N_s!}{k!l!\left(N_s-k-l\right)!}{p_1}^k{p_2}^l(1-p_1-p_2{)}^{N_s-k-l}.$$

Full text

Full counting statistics has emerged as a very powerful tool to characterize and obtain information about a quantum mechanical system by gaining knowledge on the full probability distribution function of an observable rather than just its expectation value. These concepts have been pioneered by Lesovik and Levitov Levitov et al. (1996) for transport measurements in nano-structures Qua ; a remarkable application being the demonstration of the fractional charge of quasiparticles in a fractional quantum Hall fluid de-Picciotto et al. (1998); Saminadayar et al. (1997). The concept of full counting statistics turns out to be very powerful in the context of cold atomic gases, where the observable of interest is the number of particles on a set of lattice sites, which is accessible with site-resolved quantum gas microscopes Bakr et al. (2009); Sherson et al. (2010). Especially, it has been applied to characterize quantum states of cold atomic gases in equilibrium Öttl et al. (2005); Cherng and Demler (2007); Belzig et al. (2007); Braungardt et al. (2008) as well as non-equilibrium states Lamacraft (2006); Braungardt et al. (2011); Gring et al. (2012).

Several cold atoms setups for fermionic atoms are currently equipped with a quantum gas microscope Haller et al. (2015); Parsons et al. (2016); Cheuk et al. (2016); Drewes et al. (2016); Mazurenko et al. (2017); Hilker et al. (2017); Mitra et al. (2017) and have achieved single-site and single-atom detection, which is required for the measurement of probability distributions in Fock space by accumulating histograms of particle configurations over independent measurement realizations Mazurenko et al. (2017). Notably, both the repulsive Parsons et al. (2016); Cheuk et al. (2016); Drewes et al. (2016); Mazurenko et al. (2017); Hilker et al. (2017) and the attractive Mitra et al. (2017) Hubbard model have been realized in cold atoms experiments at temperature scales that may be already relevant in the context of normal-state properties of high-$T_{c}$ cuprates. While for bosonic systems, the full counting statistics is accessible within path integral quantum Monte Carlo simulations, in the fermionic situation the notorious sign problem renders such path integral approaches inefficient. In turn, for auxiliary field QMC methods such as determinantal quantum Monte Carlo (DQMC) Ass , suitable to study such fermionic systems in certain parameter regimes, there is no direct correspondence between the computational configuration space of auxiliary fields and states in Fock space.

In this letter, we demonstrate a method to compute the full counting statistics (FCS) within the framework of determinantal quantum Monte Carlo simulations and provide first comparisons with experiments for fermionic cold atomic gases. The main observation for this ability is the fact that DQMC simulations decomposes the interacting fermionic system into an incoherent sum over density matrices for free fermions in an external potential Grover (2013). For such free fermions, the spectrum of the reduced density matrix required for the determination of the full counting statistics is related to the eigenvalues of a one-particle correlation function Cheong and Henley (2004); Peschel (2003). Furthermore, the generating function for the FCS of the particle number, magnetization, and the staggered magnetization is still quadratic in the fermionic field operators. This allows us to calculate the relevant trace over the exponentially large Hilbert space of fermionic Fock states in all particle number sectors as the determinant of a single-particle operator. We demonstrate that full counting statistics is a powerful tool for the detection of pairing and the size of the pairing wave function by an odd-even effect in the probability distribution for the particle number. Furthermore, we present the first comparison of the FCS in the repulsive Hubbard model at half filling and demonstrate that the experiments clearly observe signatures which are not captured by the Heisenberg model.

We start with the detailed description of the numerical methods for the determination of the FCS. The quantity of interest is the probability $P(N^{\uparrow}_{A},N^{\downarrow}_{A})$ that there are $N^{\uparrow}_{A}$ fermions with spin up and $N^{\downarrow}_{A}$ fermions with spin down on a subsystem $A$ with $N_{s}$ sites, see Fig. 1(a). The observable of interest is therefore given by the spin density operator $\hat{N}^{\sigma}_{A}=\sum_{i\in A}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$ with $\sigma\in\{\uparrow,\downarrow\}$ and eigenvalues $N^{\sigma}_{A}=0,\ldots,N_{s}$. Its distribution function is most conveniently derived by the generating function

[TABLE]

Here, $\rho_{A}=\text{Tr}_{B=\bar{A}}\left(\rho\right)$ denotes the reduced density matrix with $B=\bar{A}$ the complement of subsystem $A$. Then, the joint probability $P(N^{\uparrow}_{A},N^{\downarrow}_{A})$ is determined as the coefficient of the Fourier series of the generating function

[TABLE]

with $\phi^{\sigma}_{n}=\frac{2\pi n}{N_{s}+1}$ . Note that the maximal number of fermions with a given spin in the subsystem A is limited by $N_{s}$, and therefore, there are only $N_{s}+1$ coefficients in the Fourier series. In general, one is mostly interested in the probability distributions $P_{N}(N_{A})$ for the total particle number $N_{A}=N_{A}^{\uparrow}+N_{A}^{\downarrow}$ and $P_{M}(M_{A})$ for the magnetization $M_{A}=N_{A}^{\uparrow}-N_{A}^{\downarrow}$ on subsystem A. These quantities derive directly from the joint probability distribution via

[TABLE]

with $N_{A}=0,\ldots,2N_{s}$. An analogous expression follows for $P_{M}(M_{A})=\sum_{N_{A}^{\downarrow}=0}^{N_{s}}P(M_{A}+N_{A}^{\downarrow},N_{A}^{\downarrow})$ with $M_{A}=-N_{s},\ldots,N_{s}$.

In the following, we demonstrate that we can efficiently determine the generating function $\chi$ within DQMC. The standard procedure of DQMC Ass discretizes the inverse temperature $\beta$ and then decouples the interactions with a Hubbard-Stratonovich (HS) transformation at the expense of introducing an auxiliary field at every site and time slice. Then, the partition function can be written as a sum over free fermion systems coupled to an (imaginary) time-dependent Ising field. The expectation value of an observable $\hat{O}$ is

[TABLE]

where $Z=\sum_{\{{\bf s}\}}w_{\{{\bf s}\}}$ is the partition sum and $w_{\{{\bf s}\}}$ is the weight of one auxiliary field configuration $\{\bf s\}$. The latter takes the form Ass

[TABLE]

where $\mathcal{G}^{\sigma}_{\{{\bf s}\}}=\langle\hat{c}_{i,\sigma}\hat{c}_{j,\sigma}^{\dagger}\rangle_{\{{\bf s}\}}$ denotes the single-particle Green’s function for spin species $\sigma$. Such a decomposition naturally carries over to more complex quantities such as the generating function Grover (2013)

[TABLE]

with $\chi^{\sigma}_{\{{\bf s}\}}(\phi^{\sigma})=\text{Tr}\left(\rho^{\sigma}_{\{{\bf s}\},A}\,e^{i\phi^{\sigma}\hat{N}^{\sigma}_{A}}\right)$ and $\rho^{\sigma}_{\{{\bf s}\},A}$ the reduced density matrix for the free fermions with auxiliary field configuration $\{\bf s\}$ and spin $\sigma$. It is important to stress that within a Hubbard-Stratovonich configuration the reduced density matrix factorices into a spin up and spin down part. Each can be written in the form of a “Boltzmann weight” as Chung and Peschel (2001); Cheong and Henley (2004); Peschel (2003)

[TABLE]

with the entanglement Hamiltonian for spin species $\sigma$

[TABLE]

and normalization $\mathcal{K}^{\sigma}_{\{{\bf s}\},A}=\det\left(\mathbb{1}-G^{\sigma}_{\{{\bf s}\},A}\right)$. Here, $\left[G^{\sigma}_{\{{\bf s}\},A}\right]^{ij}=\langle\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{j,\sigma}\rangle_{\{{\bf s}\};i,j\in A}$ is the one-body density matrix (OBDM) for a given HS field configuration $\{{\bf s}\}$ and with sites $i$ and $j$ restricted to subsystem $A$ Chung and Peschel (2001). Therefore, the generating function for a fixed auxiliary field configuration $\{\bf s\}$ reduces to

[TABLE]

where $\{\epsilon^{\sigma}_{\alpha}\}$ are the eigenvalues of the entanglement Hamiltonian $H^{\sigma}_{A}$ and $\{\lambda^{\sigma}_{\alpha}\}$ those of the OBDM. In deriving the above expression, we have used that the particle number operator $\hat{N}^{\sigma}_{A}=\sum_{i\in A}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$ is also a quadratic operator, and commutes with the entanglement Hamiltonian $H^{\sigma}_{A}$, i.e., they have a common eigenbasis. Within this basis the grand canonical trace can be performed analytically.

Finally, the eigenvalues $\{\epsilon^{\sigma}_{\alpha}\}$ and $\{\lambda^{\sigma}_{\alpha}\}$ are related through $\epsilon^{\sigma}_{\alpha}=\log(\frac{1}{\lambda^{\sigma}_{\alpha}}-1)$ Peschel (2003), and we obtain the important result

[TABLE]

which allows for the efficient determination of the generating function $\chi$ by quantum Monte Carlo sampling of the auxiliary field configurations. If the transformations (2) and (3) are performed in every Monte-Carlo measurement step, error bars for $P_{N}(N_{A})$ and $P_{M}(M_{A})$ can be obtained in the standard way. Note that the equal-time OBDM $G^{\sigma}_{\{{\bf s}\},A}(\tau)$ depends explicitly on imaginary time $\tau$, which is suppressed in our notation. Due to the cyclic property of the trace in Eq. (4) all imaginary time slices are equivalent and it is possible to average over them to acquire additional statistics.

In the following, we demonstrate our approach for the determination of the full counting statistics in the two-dimensional Fermi-Hubbard model Micnas et al. (1990) on a square lattice. The Hamiltonian is given by

[TABLE]

The fermionic operators $\hat{c}_{{\bf r},\sigma}$ obey canonical commutation relations $\{\hat{c}^{\dagger}_{{\bf r},\sigma},\hat{c}_{{\bf r}^{\prime},\sigma^{\prime}}\}=\delta_{{\bf r},{\bf r}^{\prime}}\delta_{\sigma,\sigma^{\prime}}$, and $\hat{n}_{{\bf r},\sigma}=\hat{c}^{\dagger}_{{\bf r},\sigma}\hat{c}_{{\bf r},\sigma}$ are the density operators while $t$ denotes the nearest neighbour hopping and $U$ the on-site interation. The model can be efficiently simulated with DQMC for the attractive case $U<0$ at arbitrary fillings, and for repulsive interactions at half-filling. Especially, at half filling and on bipartite lattices, a spin-down particle-hole transformation $\hat{c}_{{\bf r},{\downarrow}}\rightarrow(-1)^{r_{x}+r_{y}}\hat{c}^{\dagger}_{{\bf r},{\downarrow}}$ maps the attractive Hubbard model into the repulsive one.

We start with the presentation of the results for the attractive Hubbard model. The full counting statistics $P_{N}(N_{A})$ for the particle number $N_{A}$ on the sublattice A is derived within the setup illustrated in Fig. 1(a). The finite-temperature phase diagram Paiva et al. (2004, 2010) of (11) features, away from half filling, a Berezinskii-Kosterlitz-Thouless (BKT) transition at temperature $T_{\text{BKT}}$ to a quasi long-range $s$-wave superconducting state. At half-filling, there is a degeneracy of s-wave superconductivity and charge density wave order, and the $SU(2)$ symmetry of the order parameter suppresses the transition temperature to zero according to the Mermin-Wagner theorm. In both cases, there is a second clearly separated temperature scale $T^{\star}\sim|U|$ which marks the onset of pair formation without long-range phase coherence. In Fig. 1(b-d), we present the behavior of the full counting statistics in the cross-over from the BCS-type superfluid of large overlapping Cooper pairs for weak interactions to a BEC of hardcore bosonic on-site pairs Eagles (1969); Leg ; Nozières and Schmitt-Rink (1985) for strong interactions. The temperature is chosen well below the characteristic temperature $T^{\star}$ for pair formation.

We observe a strong even-odd effect for increasing interactions. We will argue in the following that this phenomenon can be well understood through the size of the pairing wave function: in the extreme limit of very strong interactions, all fermions are paired up with a pair wave function localized on a single lattice site. Then, we expect a vanishing probability to find an odd number of fermions on subsystem A. An odd number of fermions can only appear due to unpaired fermions or a pairing function spreading over several lattice sites. The latter provides only a contribution for pairs along the boundary of subsystem A. In order to quantify this effect, we introduce $P_{\text{odd}}=\sum_{N_{A}\text{odd}}P(N_{A})$ as a measure of the even-odd splitting. We expect a scaling behavior

[TABLE]

with increasing subsystem size $L_{A}$ and atomic density $\langle n\rangle$; $4L_{A}\cdot\langle n\rangle$ estimates the number of pairs along the boundary. This ansatz is motivated by the idea of independent pairs randomly distributed in the sample. Then, for a pair centered next to the boundary, $p_{\text{deloc}}$ denotes the weight of the pair wave function to be on the other side of the boundary. Therefore, we can interpret $p_{\text{deloc}}$ as a measure for the size of the wave function. Indeed, we find a perfect fitting of $P_{\text{odd}}$ to (12) for a large parameter regime, see dashed lines in Fig. 2(a), which allows us to extract the value $p_{\text{deloc}}$ for different interactions and densities. The resulting $p_{\text{deloc}}$ for the Hubbard model are shown in Fig. 2(b) at density $\langle n\rangle=1$ and as a function of $|U|$; note that sufficiently low temperatures are required to suppress thermal pair breaking. In 2D the binding energy of a bound state is exponentially small in the attractive interaction Landau and Lifshitz (1987); Randeria et al. (1990). Therefore for $U/t=-3$, $L=12$, a temperature as low as $\beta t=24$ was necessary to achieve a fit to Eq. (12), whereas for large $|U|$, e.g. $U/t=-9$, “low enough” temperature is $\beta t=4$. In Fig. 2(c), the behavior of $p_{\text{deloc}}$ is shown for decreasing densities. The minimum of the pair size at half filling, $\langle n\rangle=1$, is a lattice effect (see Randeria et al. (1990) for the 2d continuum model), due to the logarithmically diverging van Hove singularity in the tight-binding density of states on the square lattice. As Fig. 2(c) demonstrates, BCS theory at T=0 gives a reasonably accurate estimate of the pair size even for intermediate values of $U$ and densities $\langle n\rangle$.

From this quantitative analysis of the even-odd splitting we conclude that the FCS can serve as a powerful tool to observe the formation of pairing in a superconducting state and provides useful information about the size of the pairing wave function. We expect that this analysis carries over to systems where DQMC simulations are hard, such as the repulsive Hubbard model away from half-filling, and can serve as a powerful experimental detection tool.

Finally, we study the FCS for the staggered magnetization in the repulsive Hubbard model at half-filling. Recently, this model was extensively studied experimentally using cold atoms in optical lattices Mazurenko et al. (2017) at interactions $U/t=7.2$ with a circular central region of homogeneous density involving approximately 80 sites surrounded by a dilute particle bath. Given the single-site resolution in the experiments Parsons et al. (2016), histograms of the staggered magnetization $M^{st}=\sum_{i}(-1)^{{\bf i}}(n_{i,\uparrow}-n_{i,\downarrow})$ inside the circular region were accumulated over more than 250 experimental realizations Mazurenko et al. (2017). Here, we show the original data points of Ref. Mazurenko et al. (2017) as the red histograms in Fig. 3. As the generating term for the staggered magnetization no longer commutes with the entanglement Hamiltonian in Eq. (8) our method to determine the FCS has to be modified: it is required to explicitly determine the entanglement Hamiltonian and subsequently diagonalize the Hamiltonian with the staggered term of the generating function added; the details of this procedure are discussed in the supplement material. For optimal comparison with the experiment, we determine the FCS of the staggered magnetization inside the same circular geometry as in the experiment, which is embedded in a large system of size $L\times L$ with $L=20$ and periodic boundary conditions. As shown in Fig. 3, we find excellent agreement without any adjustable fit parameter. Furthermore, we find strong deviation from the predictions of the staggered magnetization from the $S=1/2$ antiferromagnetic Heisenberg model. The latter takes the form $\mathcal{H}_{SO}=J\sum_{\langle i,j\rangle}\hat{{\bf S}}_{i}\cdot\hat{{\bf S}}_{j}$ with super-exchange coupling $J=4t^{2}/U$, and follows as limiting theory of the Hubbard model for large interactions $U\gg t$; the corresponding temperature is given by $T/J\,$=$\,(T/t)U/(4t)$. The Heisenberg model is simulated at the corresponding temperatures with the same circular geometry of $N_{s}=80$ sites which is embedded in a larger system of linear dimension $L=32$. For all temperatures, the half width of the distribution of the staggered magnetization for the Hubbard model is systematically smaller than that of the Heisenberg model. We expect charge fluctuations, which already play an important role at the given interaction strength $U/t=7.2$, to be the cause of this difference: it is well established, that the next order correction to the Heisenberg model in $t/U$ leads to second and third neighbour exchange as well as four-spin ring exchange interactions of order $t^{4}/U^{3}$ MacDonald et al. (1988); Takahashi (1977). Furthermore, the operator measuring the staggered magnetiztion is reduced by a renormalization factor $\left(1-8\frac{t^{2}}{U^{2}}\right)$ due to doublon-hole pairs Delannoy et al. (2005). We therefore conclude that at the intermediate interaction strength $U/t=7.2$ the differences between the Hubbard model and the Heisenberg model already play a significant role. The perfect agreement of the experimental data with our DQMC analysis demonstrates that the experiments are indeed capable of probing these corrections to the Heisenberg model.

However, there is an important comment of caution in order for the comparison with our theoretical predictions and the experimental results: As a result of the parity projection, the measurement technique of Parsons et al. (2016); Mazurenko et al. (2017) counts doublons and holes sitting on the same sublattice (i.e. on diagonally opposite corners of a plaquette) incorrectly to the staggered magnetization. On the other hand, doublon-hole pairs on neighbouring sites, which are the dominant charge fluctuations, are counted correctly to the staggered magnetization as a consequence of the staggering factor. Therefore, the experiments are currently exactly performing the measurement to observe the leading deviations between Heisenberg and Hubbard model. However, the observation of further corrections for lower interactions will require the precise experimental detection of the staggered magnetization.

In conclusion, we present a method to compute full quantum mechanical probability distributions of quadratic operators in interacting fermion systems which can be simulated with DQMC. We find that for the attractive Hubbard model the dependence of an even-odd splitting in the particle number distribution function allows one to infer the size of a preformed pair or Cooper pair from in situ images. Furthermore, we apply our method to the repulsive Hubbard model, which has recently been studied experimentally Mazurenko et al. (2017). The excellent agreement with our analysis demonstrates that the experiments are capable of observing corrections to the Heisenberg model for the intermediate interaction strengths $U/t=7.2$. Our method is also suitable for the study of the high-temperature phase of the repulsive Hubbard model away from half filling or the investigation of the universality of distribution functions Lovas et al. (2017) at a quantum critical point Meng et al. (2010); Sorella et al. (2012).

Acknowledgments

We acknowledge M. Greiner for providing the raw experimental data for Figure 3. S.H. thanks T. Roscilde and Feiming Hu for fruitful discussions as well as S. Wessel for providing benchmark results in one dimension and D. Huerga for comments on the manuscript. Computations were performed on JURECA, Jülich Supercomputing Center. Support by the German Science Foundation (DFG) through SFB TRR21 is acknowledged.

I Appendix

I.1 FCS of the staggered particle number

To compute the full counting statistics (FCS) of the staggered particle number $N_{A}^{st}=\sum_{i\in A}(-1)^{{\bf i}}(\hat{n}_{i,\uparrow}+\hat{n}_{i,\downarrow})$ (which corresponds to the staggered magnetization plus a constant, $N_{A}^{(st,U<0)}\rightarrow M_{A}^{(st,U>0)}=\sum_{i\in A}(-1)^{{\bf i}}(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})+\sum_{i\in A}(-1)^{{\bf i}}$, in the half-filled repulsive Hubbard model) a modification of the method described in the main text is necessary. The generating function of $P(N_{A}^{st})$ in one Hubbard-Stratonovich sample is

[TABLE]

Unlike the particle number $\hat{N}_{A}^{\sigma}=\sum_{i\in A}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$, the operator $i\phi\sum_{i\in A}(-1)^{{\bf i}}\hat{c}_{i,\sigma}^{\dagger}\hat{c}_{i,\sigma}$ does not commute with the single-particle entanglement Hamiltonian $H_{A}^{\sigma}$ due to the staggering factor $(-1)^{{\bf i}}$. Therefore, there is no common eigenbasis in which one can simply add the eigenvalues of the two operators as was done in the last step leading to Eq. (8) of the main text. Instead, it is necessary to compute $H_{A}^{\sigma}$ explicitly in the site basis, add the non-commuting operators

[TABLE]

and diagonalize the resulting modified entanglement Hamiltonian $\tilde{H}_{A}^{\sigma}$, which gives the eigenvalues $\tilde{\epsilon}^{\sigma}(\phi)$. The peculiarity of DQMC that the equal-time Green’s function $G^{\sigma}_{\{{\bf s}\},A}$ is non-Hermitian and does not necessarily have a spectral decomposition leads to the complication that the entanglement Hamiltonian $H_{A}^{\sigma}$ in Eq. (7) of the main text,

[TABLE]

has to be computed explicitly through the power series of the matrix logarithm Al-Mohy and Higham (2012). This needs to be done only once per Hubbard-Stratonovich sample as $H_{A}^{\sigma}$ can be reused in (14) for different values of $\phi$.

The grand-canonical trace (13) over fermionic degrees of freedom results in a single-particle determinant which is expressed in terms of the eigenvalues $\tilde{\epsilon}^{\sigma}(\phi)$ as

[TABLE]

The computation of the determinant (16) may be severely affected by numerical inaccuracies. For low temperatures the occupation numbers $\lambda^{\sigma}_{\alpha}$ in the free fermion system tend to zero or one and, consequently, the OBDM $G^{\sigma}_{\{{\bf s}\},A}$ is a nearly singular matrix. The evaluation of (15) in finite-precision arithmetic is beset with numerical instabilities and, even if a high-quality implementation of the matrix logarithm Al-Mohy and Higham (2012) is used, not all eigenvalues $\epsilon^{\sigma}_{\alpha}$ of the entanglement Hamiltonian $H_{A}^{\sigma}$ can be obtained with sufficient accuracy across the entire spectrum. Consider small (negative) “entanglement energy” $\epsilon^{\sigma}_{\alpha}$ for which according to the Fermi-Dirac statistics $\lambda^{\sigma}_{\alpha}=\frac{1}{1+e^{\epsilon^{\sigma}_{\alpha}}}\approx 1$. The evaluation of $\epsilon^{\sigma}_{\alpha}=\log\left(\frac{1}{\lambda^{\sigma}_{\alpha}}-1\right)\approx\log(1-1)$, or equivalently the evaluation of (15), is very inaccurate for the small (negative) real part of the entanglement spectrum since $\log(x)$ varies greatly for $x\rightarrow 0^{+}$. In view of round-off and cancellation errors it is therefore crucial to compute (15) with the mathematically equivalent formula

[TABLE]

This expression shifts the accuracy to small (negative) $\epsilon^{\sigma}_{\alpha}$ (large $\lambda^{\sigma}_{\alpha}\approx 1$) since $\log(x)$ is not much affected by errors in its large argument $x=\frac{\lambda^{\sigma}_{\alpha}}{1-\lambda^{\sigma}_{\alpha}}$. For the computation of (16) errors in the large (positive) part of the entanglement spectrum $\epsilon^{\sigma}_{\alpha}$ are not problematic because large $\tilde{\epsilon}^{\sigma}_{\alpha}(\phi)$ are irrelevant in the factor $(1-\lambda^{\sigma}_{\alpha})\left(1+e^{-\tilde{\epsilon}^{\sigma}_{\alpha}(\phi)}\right)$. Here it is assumed that the spectrum $\tilde{\epsilon}^{\sigma}_{\alpha}(\phi)$ of the modified entanglement Hamiltonian (14) is qualitatively similar to that of the entanglement Hamiltonian (15) in the sense that the accuracy of the relevant small $\tilde{\epsilon}^{\sigma}_{\alpha}(\phi)$ is not affected by the inaccuracy of large $\epsilon^{\sigma}_{\alpha}$.

I.2 FCS for a BCS mean-field state

The Hamiltonian of the single-band Fermi-Hubbard model in momentum space reads

[TABLE]

with the single-particle band structure $\varepsilon_{{\bf k}}=-2t(\cos(k_{x}a)+\cos(k_{y}a))$ for the square lattice. The standard BCS mean-field analysis starts with the BCS reduced Hamiltonian which neglects the sum over ${\bf q}$ in (17) and considers only scattering between pairs of zero center-of-mass momentum (${\bf q}=0$). Under the assumption of a non-zero $s$-wave pairing order parameter $\Delta=\sum_{{\bf k^{\prime}}}\langle\hat{c}_{{\bf k},\uparrow}\hat{c}_{-{\bf k},\downarrow}\rangle$ a mean-field decoupling is performed. The resulting quadratic Hamiltonian is then solved with a Bogoliubov transformation Schrieffer (1964) with coefficients

[TABLE]

which are parametrized by the chemical potential $\mu$ and $\Delta$, the gap parameter. $E_{{\bf k}}^{2}=(\varepsilon_{{\bf k}}^{(\text{HF})}-\mu(n,U))^{2}+\Delta^{2}(n,U)$ are the excitation energies of the fermionic Bogoliubov quasiparticles. In order to take into account density-density interactions at the mean-field level, it is necessary to include the Hartree-Fock potential in the single-particle energies: $\varepsilon_{{\bf k}}^{(\text{HF})}=\varepsilon_{{\bf k}}-|U|\frac{n}{2}$ where $n=(N_{\uparrow}+N_{\downarrow})/N_{\text{sites}}$ is the filling. The mean-field ground state takes the standard form in terms of the BCS coefficients (18)

[TABLE]

where $\mu(n,U)$ and $\Delta(n,U)$ are self-consistent solutions Leg of the gap equation

[TABLE]

and the number equation

[TABLE]

For general filling $n$ and Hubbard interaction $U$ these self-consistent equations need to be solved numerically to obtain $\mu$ and $\Delta$. At half filling ($n=1$) the number equation has the solution $\mu=-|U|/2$ for any choice of $\Delta$ provided that the single-particle dispersion relation $\varepsilon_{{\bf k}}$ is particle-hole symmetric around $\varepsilon_{|{\bf k}|=k_{F}}=0$ such that the integral in (21) vanishes. Thus, the inclusion of the Hartree shift in the single-particle energies has preserved, at the mean-field level, the particle-hole symmetry which the Hubbard model posseses for half filling. In the atomic limit $t/U=0$ the solutions are $\mu=-|U|/2$ and $\Delta=|U|\sqrt{n(2-n)}/2$.

The BCS state is the ground state of a quadratic (mean-field) Hamiltonian and therefore Wick’s theorem can be used to factorize any correlation function into sums of products of single-particle Green’s functions. Following closely the method of Ref. Cherng and Demler (2007), the generating function for the total particle number is written as

[TABLE]

where $E_{i\sigma}=\hat{c}_{i\sigma}+\hat{c}_{i\sigma}^{\dagger}$ and $F_{i\sigma}=\hat{c}_{i\sigma}+e^{i\phi}\hat{c}_{i\sigma}^{\dagger}$. In the second equation of (22) the particle-hole symmetry of the Hubbard Hamiltonian was used to replace the total particle number by the total hole number, which simplifies the analysis Cherng and Demler (2007) as the identity $e^{i\lambda(1-\hat{c}_{i}^{\dagger}\hat{c}_{i})}=E_{i}F_{i}(\phi)$ can then be used. The multipoint correlation function (22) is contracted according to Wick’s theorem, and, taking into account the fermionic anticommutation relations, the non-zero full contractions Cherng and Demler (2007) constitute a determinant:

[TABLE]

In the last step we have introduced the eigenvalues $\mu_{k}$ of the matrix of normal and anomalous Green’s functions:

[TABLE]

The Green’s functions, which are restricted to subsystem $A$, are computed for the BCS state (19). As described in detail in Cherng and Demler (2007), complex conjugate pairs of eigenvalues of $M$ can interfere in Eq. (24) leading to a suppression of odd versus even values in the particle number distribution $P(N_{A})$. Complex eigenvalues can appear as soon as non-zero values of the anomalous Green’s functions $\langle\hat{c}_{i\uparrow}\hat{c}_{j\downarrow}\rangle_{i,j\in A}=-\langle\hat{c}_{j\downarrow}\hat{c}_{i\uparrow}\rangle_{i,j\in A}\neq 0$ destroy the hermiticity of $M$, which makes evident that BCS pairing is at the origin of the even-odd effect in the distibution $P(N_{A})$.

With increasing $|U|$ the local number fluctuations of the BCS state (19) acquire a large unphysical extensive (Poissonian) contribution (see Figs. 1(c) and (e) of the main text), which is partly due to the fact that the mean-field Hamiltonian does not conserve the total number of particles. Furthermore, the wave function (19), when applied to a lattice model, fails to account for the nearest-neighbour repulsion Micnas et al. (1990) of tightly bound pairs in the limit of strong correlations $|U|\gg t$, where interactions between Bogoliubov quasiparticles need to be included Belkhir and Randeria (1992). This is different from the equivalent model in the continuum where (19) becomes the exact ground state wavefunction both in the extreme BCS and BEC limit Eagles (1969); Leg ; Nozières and Schmitt-Rink (1985) and reproduces the correct particle number variance and higher cumulants Belzig et al. (2007).

I.3 Trinomial distribution

A simple model that is capable of describing the distribution of the staggered magnetization $M_{A}^{st}$ in the limit of high temperatures (see Fig. 3(d) of the main text) regards the staggered magnetization at each lattice site as an independent random variable $m_{i}^{st}$ which can take on the values $0,+1,-1$ with the probabilities

[TABLE]

with $p(\sigma,f)$ denoting the probability of the elementary event that a spin $\sigma\in(\uparrow,\downarrow)$ is placed on a lattice site with staggering factor $(-1)^{{\bf i}}\equiv f\in(+,-)$ and $d$ and $h$ signifying the placement of a doublon and hole, respectively. In the presence of particle-hole symmetry and equal chemical potentials for spin up and down, $p(d)=p(h)$ and $p(\uparrow,f)=p(\downarrow,-f)$, such that the probabilities of all three elementary events are described by the single parameter $p_{d}$, which we set to the average double occupancy $p_{d}\leftarrow\langle d\rangle=\frac{1}{N}\sum_{i}\langle\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}\rangle$ as computed with Monte Carlo. The total staggered magnetization on a subsystem with $N_{s}$ sites is given by a sum over a trinomial distribution

[TABLE]

where

[TABLE]

For the distribution function in Fig. 3(d) of the main text, the value $\langle d\rangle=0.066(1)$ as extracted from the Monte Carlo simulations at temperature $T/t=1.14$ was used.

I.4 Additional datasets: FCS of magnetization, small subsystems

The FCS of the magnetization $\hat{M}_{A}=\sum_{i\in A}(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})$ (see Fig. A4) agrees much better between Hubbard and Heisenberg model than that of the staggered magnetization. This can be traced back to the fact that the total magnetization $\hat{M}^{\text{tot}}=\sum_{i}(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow})$ commutes with the Hubbard Hamiltonian $\mathcal{H}$ and the magnetization on a subsystem $A$ commutes up to an operator $\hat{O}_{\partial A}$ that has support only on the boundary of $A$, $[\mathcal{H},\hat{M}_{A}]=\hat{O}_{\partial A}$, which comes from the change of magnetization due to particles hopping into and out of the subsystem. A possible renormalization factor of the magnetization that derives from the canonical transformation of the Hubbard model into a spin-only Hamiltonian Delannoy et al. (2005) should then be at most a boundary effect. For temperatures $T/t=0.25$ and $T/t=0.35$ a tiny even-odd modulation is visible in the FCS of the magnetization. This is the precursor of the full suppression of odd magnetizations in the FCS of the Heisenberg model. Unlike in Fig. 3 of the main text, in Fig. A4 odd magnetizations have not been binned together with the nearest even magnetization to give a smooth histogram. Instead, the probability distribution $P(M_{A})$ for the Heisenberg model has been divided by the bin width, i.e. by two, to aid the visual comparison with the probability distribution for the Hubbard model which has twice as many possible values for $M_{A}$ on the abscissa.

In Fig. A5 the distribution $P(M_{A})$ of the magnetization in the repulsive Hubbard model at half filling is shown for small subsystems of size $3\times 3$ and $4\times 4$. The interaction strength ranges from $U=4$ to $U=14$ and the inverse temperature from $\beta t=0.5$ to $\beta t=2$. Already at these high temperatures an even-odd splitting is visible for large enough Hubbard repulsion $U$.

Fig. A6 shows the distribution $P(M^{st}_{A})$ of the staggered magnetization $\hat{M}^{st}_{A}=\sum_{i\in A}(-1)^{{\bf i}}(\hat{n}_{i,\uparrow}-\hat{n}_{i,\downarrow})$ in the half-filled repulsive Hubbard model at $U/t=7.2$ for small subsystems $A$. In the left column of Fig. A6, the FCS is shown for the highest and lowest temperatures realized in the experiment of Ref. Mazurenko et al. (2017). For the highest temperature $T/t=1.14$, we find good agreement with a simple model based on the trinomial distribution (magenta line) which treats each lattice site as independent (see previous section) and is parametrized solely by the average double occupancy $p_{d}=\frac{1}{N}\sum_{{\bf r}}\langle\hat{n}_{{\bf r},\uparrow}\hat{n}_{{\bf r},\downarrow}\rangle$ = 0.066(1) as computed with Monte Carlo for the given temperature. Interestingly, the even-odd splitting in $P(M^{st}_{A})$ is smeared out much more quickly with increasing subsystem size than in the case of $P(M_{A})$. The origin of the even-odd asymmetry is in both cases the formation of local moments; however, contrary to $M_{A}$ ($N_{A}$ for $\mathcal{H}_{U<0}$ and half filling), which can only change by $\Delta M_{A}=\pm 1$ when a particle crosses the boundary of subsystem $A$, $M^{st}_{A}$ ($N^{st}_{A}$ for $\mathcal{H}_{U<0}$ and half filling) is additionally affected by hopping processes in the bulk of $A$, which change the staggered magnetization by $\Delta M^{st}_{A}=\pm 2$. At this level of analysis and at half filling, the even-odd splitting in $P(M^{st}_{A})$ provides no additional information beyond the parameter $p_{\text{deloc}}$ extracted earlier from the even-odd splitting in $P(M_{A})$ (or rather $P(N_{A})$ for $H_{U<0}$).

I.5 Benchmarking

In order to verify the correctness of our numerical implementation we compare our DQMC method with Stochastic Series Expansion Sandvik (1999) (SSE) QMC for a Fermi-Hubbard chain Sandvik (1992), where the FCS can be obtained simply by accumulating histograms over Monte Carlo measurement steps. Open boundary conditions are necessary to ensure that the SSE QMC has no sign problem. The agreement between the two methods is excellent (see Fig. A7).