Critical Opalescence across the Doping Driven Mott Transition in Optical Lattices of Ultracold Atoms

TL;DR

This paper investigates critical opalescence in ultracold atom optical lattices, revealing fluctuation behavior near a phase transition in the fermionic Hubbard model, with implications for understanding high-temperature superconductivity.

Contribution

It provides a detailed theoretical analysis of density fluctuations across the Mott transition, highlighting the role of nonlocal correlations and predicting observable signatures in ultracold atom experiments.

Findings

Density fluctuations increase near the phase transition at low temperatures.

Nonlocal contributions to fluctuations decrease unexpectedly as criticality approaches.

Theoretical results align with experimental data in the high-temperature regime.

Abstract

Phase transitions and their associated crossovers are imprinted in the behavior of fluctuations. Motivated by recent experiments on ultracold atoms in optical lattices, we compute the thermodynamic density fluctuations $\delta N^2$ of the two-dimensional fermionic Hubbard model with plaquette cellular dynamical mean-field theory. To understand the length scale of these fluctuations, we separate the local from the nonlocal contributions to $\delta N^2$. We determine the effects of particle statistics, interaction strength $U$, temperature $T$ and density $n$. At high temperature, our theoretical framework reproduces the experimental observations in the doping-driven crossover regime between metal and Mott insulator. At low temperature, there is an increase of thermodynamic density fluctuations, analog to critical opalescence, accompanied by a surprising reduction of the absolute value of their nonlocal contributions. This is a precursory sign of an underlying phase transition between a pseudogap phase and a metallic phase in doped Mott insulators, which should play an important role in the cuprate high-temperature superconductors. Predictions for ultracold atom experiments are made.

Figures

Equations

$$H=-\sum _{⟨ij⟩\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }-\mu \sum _{i\sigma }{n}_{i\sigma },$$

$$\delta N^2=A{k}_{B}T{n}^2{\kappa }_T$$

$$k_{B}T\frac{\partial ⟨n(\mathbf{x})⟩}{\partial \mu }=\int d{\mathbf{x}}^{\prime }\left[⟨n(\mathbf{x})n({\mathbf{x}}^{\prime })⟩-⟨n(\mathbf{x})⟩⟨n({\mathbf{x}}^{\prime })⟩\right].$$

$$n^2{\kappa }_T$$

$$=\frac{1}{A{k}_{B}T}\left[\delta n^2+\sum _{i\ne j}\left(⟨n_i{n}_j⟩-⟨n_i⟩⟨n_j⟩\right)\right],$$

$$\delta N^2=A{k}_{B}T{n}^2{\kappa }_T=\delta n^2+\delta n_{nl}^2,$$

$$\delta n^2\equiv ⟨n_i^2⟩-⟨n_i{⟩}^2=n+2D-n^2,$$

$$\delta n_{nl}^2=\delta N^2-\delta n^2=A{k}_{B}T{n}^2{\kappa }_T-\delta n^2$$

Full text

I Introduction

Theoretical models that can give, with minimal assumptions, insights on how complex behaviors and structures emerge in natural phenomena are particularly appealing. The Hubbard model for strongly interacting electrons on a lattice is such a model. It contains the competition between a kinetic energy term, that describes the delocalisation of particles on a lattice, and an on-site repulsive interaction $U$ term, that favors instead particle localisation. As a prominent example, this model has long helped to describe the Mott transition Mott (1974); Imada et al. (1998), where electron-electron repulsion turns a metal into an insulator. It has also helped to throw light on how complex phases of matter in quantum materials, such as cuprate high-temperature superconductors Anderson (1987); Lee et al. (2006); Tremblay et al. (2006), arise from strong interactions.

Nevertheless, correlated electron systems pose fundamental challenges to many-particle quantum theory since even a minimal model such as the Hubbard model is hard to solve or simulate Troyer and Wiese (2005). Ultracold atoms in optical lattices have emerged as a useful platform for implementing the Hubbard model and thereby for performing “quantum simulations” of interacting quantum systems Gross and Bloch (2017); Bloch et al. (2008); Jaksch and Zoller (2005); Lewenstein et al. (2007); Esslinger (2010). The synergy between advances such as single atom microscopy Bakr et al. (2009); Sherson et al. (2010); Cheuk et al. (2015); Haller et al. (2015); Parsons et al. (2015); Omran et al. (2015); Edge et al. (2015) in ultracold atom experiments and theoretical research is leading to a better understanding of the phase diagram of the Hubbard model. Important milestones include the observation of the Mott insulating phase Jordens et al. (2008); Schneider et al. (2008); Cheuk et al. (2016a); Greif et al. (2016), the measure of the equation of state Cocchi et al. (2016); Hofrichter et al. (2016), the detection of short-range spin and charge correlations Cheuk et al. (2016b); Parsons et al. (2016); Boll et al. (2016); Drewes et al. (2016, 2017) and their implication for transport properties Nichols et al. (2018); Brown et al. (2019), the measure of entanglement Islam et al. (2015); Cocchi et al. (2017); Walsh et al. (2019a, b), and the observation of long-range antiferromagnetism Mazurenko et al. (2017).

A central programme of the synergic effort in theory and experiments with ultracold atoms is the investigation of the Hubbard model in two dimensions. This is motivated by the link to the physics of cuprate high-temperature supercondutors. These materials are doped Mott insulator whose important properties should be captured by this model.

In cold atoms investigations, density fluctuations measured through single-atom microscopy have emerged as a powerful tool Müller et al. (2010); Sanner et al. (2010); Drewes et al. (2016) to detect and characterize phase transitions and crossovers, and hence to throw new light on the phase diagram of the Hubbard model. For example, experiments in Ref. Drewes et al., 2016 compared thermodynamic and local density fluctuations, unveiling a clear metal to Mott insulator crossover as a function of particle density (or doping): in the Mott insulating state thermodynamic and local fluctuations are both suppressed, whereas in the metallic state density fluctuations of all length scale, and not only local ones, contribute to the thermodynamic fluctuations.

What happens to this metal to Mott insulator crossover at much lower temperatures? While such temperatures are currently inaccessible experimentally, this may change soon. Similarly, what happens to crossovers as a function of doping? The answer to this question bears key implications for the physics of the cuprates. Upon doping the Mott insulator, a pseudogap phase emerges. The nature of this phase is under intense scrutiny in cuprates. Based on recent theoretical work Sordi et al. (2010, 2011, 2012a) with the cellular extension Maier et al. (2005); Kotliar et al. (2006); Tremblay et al. (2006) of dynamical mean-field theory Georges et al. (1996) (CDMFT), the key qualitative feature of the temperature-doping phase diagram of the two-dimensional Hubbard model is a finite-temperature and finite-doping first-order transition between a metal at low density and a pseudogap at high density Sordi et al. (2010, 2011), which ends in a second-order critical endpoint from which a sharp crossover, called Widom line, emerges up to high temperatures Sordi et al. (2012a). The defining feature of the critical endpoint is the divergence of the correlation length $\xi$, while the Widom line is the locus of the maxima of $\xi$ in the supercritical region. Now, density fluctuations are sensitive to $\xi$, and hence they are expected to be sensitive signatures of both the endpoint and of the associated Widom line.

These signatures are the focus of this paper. We calculate separately the local and nonlocal contributions to density fluctuations for the square lattice near-neighbor two-dimensional Hubbard model using CDMFT, spanning two decades in temperatures, ranging from temperatures accessible to experiment down to the lower temperatures needed to detect the first-order transition. In Section III we first analyze the density fluctuations at temperatures accessible to experiments, and then compare them with the available experimental data of Ref. Drewes et al., 2016. The excellent agreement between theory and experiment proves the accuracy of our approach and provides a strong framework to interpret current experiments. In Section IV, we then compute the density fluctuations at low temperatures, making predictions in Sec. V that may allow experiments to detect the Widom line, critical endpoint and first-order transition. The low-temperature predictions are robust but qualitative since a highly frustrated lattice should be used to avoid magnetic phase transitions. The fermion sign problem prevents us from computing on such lattices.

II Model and method

We consider the one-band two-dimensional Hubbard model on a square lattice. Its Hamiltonian reads

[TABLE]

where $t_{ij}$ is the hopping between first neighbors (and $t_{ij}=1$ sets our energy units), $U$ is the repulsive on-site interaction, and $\mu$ is the chemical potential. The operators $c^{\dagger}_{i\sigma},c_{i\sigma}$ respectively create and destroy a fermion with spin $\sigma=\{\uparrow,\downarrow\}$ at site $i$, and $n_{i}=c_{i}^{\dagger}c_{i}$ is the number operator. We solve this model using plaquette cellular DMFT Maier et al. (2005); Kotliar et al. (2006); Tremblay et al. (2006). This method selects a cluster of lattice sites - here we consider a $2\times 2$ plaquette - and replaces the missing lattice environment with a self-consistent bath of noninteracting fermions. We solve the cluster in a bath problem with continuous-time quantum Monte Carlo in the hybridisation expansion (CT-HYB) Gull et al. (2011). Here we focus on the normal state only, where there is an equal number of fermions with spins up and down and the spin alignments have short range correlations but no long range order. We map out the phase diagram of Eq. 1 as a function of occupation $n$ and interaction strength $U$ for two decades in temperatures, from $T\approx 2$ down to $T\approx 1/100$. 111Our CDMFT calculations required a large computation investment. At high temperature (e.g. $T=1/1.5$), one CDMFT iteration took about $20$ minutes on $48$ processors. At low temperature (e.g. $T=1/100$), the computational effort increases, and one CDMFT iteration took about 1 hour on $72$ processors. We normally reached convergence within $50$ CDMFT iterations, but close to phase boundaries and the Widom line hundreds of iterations were necessary. Once convergence is reached, averages were taken over at least the last $30$ CDMFT iterations (and up to hundreds close to phase boundaries). The resulting root mean square deviation on local quantities such as $n$ and $D$ that were needed to compute the density fluctuations were on the fifth digit.

III Density fluctuations: definitions and comparison with available experiments

To extract the density fluctuations, we follow a protocol similar to the experimental work of Drewes et al. Drewes et al. (2016). We first summarize this protocol by considering temperatures accessible to current experiments and by comparing our results to the experimental findings of Ref. Drewes et al. (2016). Then in Sec. IV that follows, we extend our analysis to lower temperatures to discuss enhanced density fluctuations in proximity to a critical point at finite $T$ and $n$ which may be detected in future experiments.

III.1 Thermodynamic fluctuations

First, consider thermodynamic density fluctuations. The simplest measure of these fluctuations is the variance of the number of fermions $N$ in the area $A$: $\delta N^{2}\equiv\langle(N-\langle N\rangle)^{2}\rangle$. The fluctuation-dissipation theorem

[TABLE]

links the thermodynamic fluctuations $\delta N^{2}$ to the charge compressibility $\kappa_{T}=\frac{1}{n^{2}}\left(\frac{\partial n}{\partial\mu}\right)_{T}$, where $n=N/A$ is the particle density. From the isothermal occupation curve $n(\mu)_{T}$ we extract $\kappa_{T}$ by numerical derivative, and hence using Eq. 2 we obtain $\delta N^{2}$.

Figure 1(a) shows the thermodynamic fluctuations $\delta N^{2}$ as a function of $n$ for different values of the interaction strength $U$. Our data are taken at temperature $T=1/1.5$, which is close to the lowest temperature that can be reached in current experiments.

We point out three interesting features. First, there is excellent agreement between our calculations (open symbols) and the experimental results of Ref. Drewes et al. (2016) (full symbols). Second, the variance $\delta N^{2}$ is much smaller than expected from the density $n$ (i.e., $\delta N^{2}/\langle N\rangle\ll 1$), implying sub-Poissonian statistics. This occurs even in the noninteracting case $U=0$. Indeed, as noticed in previous work Müller et al. (2010); Sanner et al. (2010); Drewes et al. (2016); Cheuk et al. (2016b), this suppression of the thermodynamic density fluctuations at $U=0$ is a consequence of the Pauli exclusion principle. According to this principle, two fermions cannot occupy the same quantum state, which in real space means that fermions ‘avoid each other’, or ‘antibunch’, thereby suppressing the density fluctuations. This sharply differs from the ideal classical gas, for which $\kappa_{T}=1/(nk_{B}T)$, and hence, from the fluctuation-dissipation theorem Eq. 2, $\delta N^{2}/\langle N\rangle=1$, i.e. Poissonian statistics. Third, the suppression of thermodynamic density fluctuations $\delta N^{2}$ due to Pauli principle is further enhanced by interactions. Indeed, the interaction strength $U$ further reduces $\delta N^{2}$. As expected, the closer the occupation $n$ is to half filling $n=1$, the larger the effect of $U$. Above the critical interaction $U_{\rm MIT}\sim 5.9$ where the system becomes an insulator at half-filling, the thermodynamic density fluctuations are strongly suppressed until they vanish completely at zero temperature, reflecting the incompressible nature of the Mott insulating state ($\kappa_{T}=0$).

III.2 Local and nonlocal density fluctuations

Next, we consider the length scales of the density fluctuations. From the definition of the occupation $n=n(\bf{x})$ in a translationally invariant system, one has Zhou and Ho (2011):

[TABLE]

Therefore, using the fluctuation-dissipation theorem Eq. 2, the density fluctuations at all length scales are included in $\delta N^{2}$. Eq. 3 thus suggests to write the thermodynamic fluctuations $\delta N^{2}$ as the sum of local fluctuations $\delta n^{2}$ and nonlocal fluctuations $\delta n_{nl}^{2}$ Duchon et al. (2012); Drewes et al. (2016): $\delta N^{2}=\delta n^{2}+\delta n_{nl}^{2}$. For our two-dimensional lattice, Eq. 3 reads:

[TABLE]

so

[TABLE]

where we have identified $\delta n_{nl}^{2}$ with $\sum_{i\neq j}\left(\langle n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle\right)$. Now the local (onsite) density fluctuations $\delta n^{2}$ can be calculated explicitly:

[TABLE]

where $n=n_{i}=\langle n_{i\uparrow}\rangle+\langle n_{i\downarrow}\rangle$ is the density, $\langle n_{i\uparrow}\rangle=\langle n_{i\downarrow}\rangle$ is the spin-resolved density, and $D=\langle n_{i\uparrow}n_{i\downarrow}\rangle$ is the double occupancy. So, the nonlocal fluctuations

[TABLE]

can be obtained as the difference between the thermodynamic and local fluctuations.

Figure 1(b) and Figure 1(c) show the local and nonlocal fluctuations $\delta n^{2}$ and $\delta n_{nl}^{2}$ respectively as a function of $n$ for the same model parameters as Fig. 1(a). Again our calculations (open symbols) and the experimental data of Ref. Drewes et al. (2016) (full symbols) are in excellent agreement, as shown in Figure 1(c). The behavior of the local (onsite) fluctuations $\delta n^{2}$ in Figure 1(b) can be understood by considering the $U=0$ and the $U=\infty$ limits. For $U=0$, $D=n^{2}/4$, so $\delta n^{2}=n-n^{2}/2$ (red dotted line), which is a monotonically increasing function of $n$. For $U=\infty$, $D=0$, so $\delta n^{2}=n-n^{2}$ (gray dashed line), which is a non-monotonic function of $n$. Therefore, the effect of $U$ is to strongly reduce the local density fluctuations, especially close to half filling.

By comparing the overall magnitude of $\delta N^{2}$ in Figure 1(a) and $\delta n^{2}$ in Figure 1(b), we find that the inequality $\delta N^{2}\ll\delta n^{2}$ is satisfied over the whole range of fillings $n$: clearly, the local fluctuations $\delta n^{2}$ do not account for all thermodynamic fluctuations, which also comprise nonlocal fluctuations $\delta n_{nl}^{2}$ that are negative and of comparable magnitude. Fig. 1(c) shows the filling dependence of $\delta n_{nl}^{2}$ calculated with Eq. 7. These negative contributions to $\delta N^{2}$ lead to sub-Poissonian statistics. As described in previous work Cheuk et al. (2016b); Cheuk (2017), the behavior of $\delta n_{nl}^{2}$ is controlled by the competing effects of Fermi statistics, which dominates for low density and small values of $U$, and superexchange, which conversely dominates for large values of $U$ and upon approaching half filling.

More specifically, for $U=0$, nonlocal correlations arise from the Pauli principle that keeps apart fermions with the same spin, so fermions move in a way to ‘avoid each other’, or antibunch. This causes a suppression in charge density in the vicinity of any fermion, and this charge suppression can be thought of as a hole, and is usually called exchange or Pauli hole. At low density, the interaction strength $U$ can enhance this antibunching effect, giving rise to what is known as a correlation hole Mahan (2000). As a result, from low to moderate density, $n<0.5$, $\delta n_{nl}^{2}$ grows in absolute value, i.e. becomes more negative, with increasing density since there are more fermions in the lattice and thus the average fermion spacing decreases with increasing $n$ and the antibunching effect becomes more pronounced. In that range of filling, $\delta n_{nl}^{2}$ does not vary much with increasing $U$.

By contrast, from moderate density up to half filling, $n>0.5$, $U$ has an appreciable effect. We see in Figure 1(c) that as $U$ increases, $\delta n_{nl}^{2}$ decreases in absolute value, i.e. become less negative. This is because upon approaching half filling $n=1$, the interaction strength $U$ reduces the double occupancies, so the fermions become more localised around their sites, meaning that the charge becomes frozen. In that region near half-filling, and at large $U$, another way to understand negative non-local correlations is to realize that they contribute to lower the energy through superexchange $J$. Indeed, $J$ comes from virtual double occupation produced by a fermion leaving a hole behind and hopping on a site occupied by a fermion of opposite spin. Effectively, this creates correlations between empty sites (holons) and doubly occupied sites (doublons) on neighboring sites 222Virtual processes occur across a bond between neighboring sites of the lattice. It is intriguing that the virtual process description becomes consistent near $n=0.5$, the classical bond-percolation threshold of the two-dimensional square lattice.. Holons and doublons are then close together, or ‘bunched’, i.e. they form quasi-bound states that lead to negative near-neighbor correlations as can be seen from the definition of nonlocal fluctuations: when there is a virtual double occupancy, then a neighboring site will be empty and the product $n_{i}n_{j}$ will be zero so that the contribution of that process to $\sum_{i\neq j}\left(\langle n_{i}n_{j}\rangle-\langle n_{i}\rangle\langle n_{j}\rangle\right)$ will be negative. However, since $J=4t^{2}/U$ is inversely proportional to $U$, the larger the value of $U$, the weaker this doublon-holon correlation since there is no more energy gain associated with virtual hops. Eventually then, non-local fluctuations decrease with increasing $U$. In the Mott insulating state at large $U$ and $n=1$, thermodynamic fluctuations essentially become local fluctuations. This description is complementary to the explanation in Refs. Cheuk et al. (2016b); Cheuk (2017) and Ref. Drewes et al. (2016). The doublon-holon correlations are also important in spin transport Nichols et al. (2018).

III.3 Temperature dependence of nonlocal density fluctuations

Up to now we considered one temperature only, $T=1/1.5$, which is close to the lowest temperature reached by ultracold atom experiments in Ref. Drewes et al. (2016). Further insight can be gained from the temperature evolution of the nonlocal fluctuations $\delta n_{nl}^{2}$. Figure 2 shows $\delta n_{nl}^{2}$ as a function of temperature for different values of $U$, at half-filling $n=1$. An important finding of Ref. Drewes et al. (2016) is that nonlocal fluctuations decrease in absolute value upon increasing $T$ (see experimental data in our Fig. 2). Physically, this is because with increasing $T$ the De Broglie thermal wavelength $\xi_{th}$ becomes smaller than the average distance between fermions, freezing nonlocal quantum fluctuations, hence recovering classical Poissonian statistics.

Our results in the temperature range $1/1.5<T<2$, displayed in Fig. 2, confirm this behavior and are in good quantitative agreement with experimental data (solid symbols). Furthermore, in our CDMFT simulations we can access lower temperatures than those experimentally available. In a metallic state, for $U<U_{\rm MIT}$, $\delta n_{nl}^{2}$ remains a monotonic function of $T$ at all temperatures and thus becomes more negative with decreasing $T$. In sharp contrast, for $U>U_{\rm MIT}$, a new temperature scale emerges where $\delta n_{nl}^{2}$ becomes less negative or saturates at low temperature. This occurs because the system is entering the Mott insulating state as the temperature is reduced. Hence, states become more localized, or from another point of view, the doublon-holon correlations lead to a reduction of the magnitude of the nonlocal correlations.

Far away from half filling, $n\leq 0.75$, the system is metallic and Ref. Drewes et al. (2016) demonstrates that $\delta n_{nl}^{2}$ is a monotonic function of $T$, as for the case $U<U_{\rm MIT}$ at $n=1$ shown in Fig. 2. Thus in a metallic state it is the De Broglie thermal wavelength $\xi_{th}$ that controls the temperature behavior of $\delta n_{nl}^{2}$. In the next section we show that, for $U>U_{\rm MIT}$ and $n>0.75$ (so in a doped Mott insulator), a novel emergent phenomenon due to the interplay of Mott physics and short-range correlations is controlling the low temperature behavior of $\delta n_{nl}^{2}$.

IV Density fluctuations across the doping driven Mott transition

In this section we analyze the behavior of the density fluctuations at temperatures lower than those reached by current experiments.

IV.1 First-order transition, critical point, and Widom line

According to previous theoretical work with plaquette CDMFT Sordi et al. (2010, 2011), the key qualitative feature of the $U-n$ phase diagram of the normal state two-dimensional Hubbard model is a low-temperature first-order transition extending from half filling $n=1$ to away from half filling, as sketched with a blue line in the $U-n$ phase diagram in Figure 3(a). At $n=1$ this transition separates a metal from a Mott insulator as $U$ increases (this corresponds to the familiar $U$ driven Mott transition). For $n\neq 1$ this transition separates two metallic states: as $n$ is lowered, one moves from a strongly correlated metal with a pseudogap to a correlated metal.

Adding the temperature axis, the first-order transition line in the $U-n$ phase diagram of Fig. 3(a) becomes a surface of first-order transition in the $U-n-T$ space. Upon increasing $T$, this surface ends in a second order critical line, from which a crossover emerges. A slice of the $U-T-n$ space at constant $U$, with $U>U_{\rm MIT}$, thus gives rise to the $T-n$ phase diagram of Fig. 3(b), which contains our CDMFT results for $U=6.2$ (i.e. a value of the interaction $U$ slightly larger than the critical threshold for obtaining a Mott insulator at $n=1$). The first-order transition is delimited by the spinodal lines $n_{c1}$ and $n_{c2}$ (line with solid triangles), where the pseudogap phase and the metallic phase cease to exist, respectively. This first-order transition ends in a second order critical endpoint at $(n_{c},T_{c})\sim(0.95,0.015)$ ($T_{c}\sim 0.004$ in units of the bare half-bandwidth), from which a rapid crossover, known as Widom line Xu et al. (2005); McMillan and Stanley (2010), emerges. The endpoint is defined as the location of the divergence of the correlation length and the Widom line as the locus of maxima of the correlation length. As $U$ increases, the critical temperature $T_{c}$ decreases rapidly while the corresponding critical filling $n_{c}$ increases. While the first-order transition may be masked by long-range ordered phases at low temperature, one expects that in frustrated systems the Widom line will continue to describe high-temperature crossovers.

Operationally we construct the Widom line with the local maxima of the isothermal charge compressibility $\kappa_{T}$ as a function of $n$ (blue dashed line in Fig. 3(b)). As shown in Fig. 3(c), the rather featureless $\kappa_{T}(n)$ at high temperature develops a broad peak that narrows and grows rapidly in intensity as $T$ decreases towards $T_{c}$, where $\kappa_{T}$ diverges. Hence, a region of enhanced charge compressibility emanates from the endpoint into the supercritical region $T>T_{c}$ (see Fig. 3(b) and Fig. 3(c)), all the way up to $T\approx 5T_{c}$ Sordi et al. (2012a). Crossing the Widom line along isotherms or along paths at constant $n$ produces sharp but continuous changes in the thermodynamic Sordi et al. (2012a) properties, such as compressibility, dynamic Sordi et al. (2012a) properties such as the local density of states, and transport properties such as $c$-axis resistivity Sordi et al. (2013). In the following, we focus on the changes of the density fluctuations upon crossing the Widom line.

IV.2 Critical opalescence

A defining feature of the critical point is the divergence of the correlation length $\xi$, which quantifies the characteristic distance over which fermion-density fluctuations are correlated. The divergence of $\xi$ is linked to the development of fluctuations in the density at all scale up to infinity, leading to the well known phenomenon of critical opalescence: physically, in our case, this means that any two fermion-density fluctuations are correlated, no matter the distance between them. Thus the thermodynamic fluctuations $\delta N^{2}$ are expected to diverge. The fluctuation-dissipation theorem associates the divergence of $\delta N^{2}$ with the divergence of $\kappa_{T}$. This expectation is fully confirmed by the calculation of $\delta N^{2}(n)$ for different temperatures at $U=6.2$, as shown in Figs. 4(a)(d). Away from half-filling, in the region $0.92<n<1$ [shaded region in Fig. 4(a)], the thermodynamic fluctuations $\delta N^{2}$ show an anomalous peak that narrows in size and increases in magnitude as the endpoint $T_{c}$ is approached from above (see Fig. 4(d)). It extends up to $T\approx 1/14$ or $T\approx 5T_{c}$. The position of these peaks are the same of those defining the Widom line. Detection of this non monotonic behavior is a precursor of the endpoint. This behavior is not surprising since we used the fluctuation-dissipation theorem Eq. 2 to calculate $\delta N^{2}$, and $\kappa_{T}$ diverges at $T_{c}$ (Fig. 3(b) and Ref. Sordi et al. (2012a) for a detailed discussion).

IV.3 Local and nonlocal fluctuations near criticality

Since the pioneering work of Kadanoff Kadanoff (1966) and Wilson Wilson (1971, 1975), it is well known that at the critical point there are fluctuations at all length scales: at $T_{c}$ the scale of the largest density fluctuations becomes infinite, however fluctuations at smaller scales persist, as illustrated by the phenomenon of critical opalescence in fluids. What remains interesting then, is to pin down the contribution of different length scales of the thermodynamic density fluctuations. Here we focus on the relative contributions of local and of nonlocal fluctuations to the thermodynamic fluctuations $\delta N^{2}$.

Let us begin with the local fluctuations $\delta n^{2}$. Figure 4(b) shows $\delta n^{2}$ versus $n$ for $U=6.2$ near $T_{c}$. In sharp contrast with $\delta N^{2}$, $\delta n^{2}(n)$ does not show an enhancement along the Widom line (see Fig. 4(e)). Its temperature dependence is controlled by that of double occupancy so it is very weak when $T$ is much less than $U$. Local fluctuations $\delta n^{2}$ are mostly dominated by quantum effects, as we discuss momentarily.

Next we turn to the nonlocal fluctuation $\delta n_{nl}^{2}$, shown in Figure 4(c) again for a temperature close to $T_{c}$. Upon crossing the Widom line, $\delta n_{nl}^{2}$ has a non-monotonic behavior, showing an additional upward peak as a function of $n$, which sharpens and narrows as the critical endpoint is approached (see Fig. 4(f)). Hence, quite surprisingly, $\delta n_{nl}^{2}$ becomes less negative along the Widom line. Note that away from the Widom line (say for $n<0.9$), $\delta n_{nl}^{2}$ becomes more negative with decreasing $T$, in agreement with the discussion of Figure 1(c).

To understand the behavior of $\delta n_{nl}^{2}$ and of $\delta n^{2}$ one needs to realize that close to the endpoint and to the Widom line, the system enters into the so-called renormalized classical regime Chakravarty et al. (1989). In this regime, the fluctuations are large and the correlation length $\xi$ becomes larger than the thermal De Broglie wavelength $\xi_{th}$, and hence the most important fluctuations are classical. A large correlation length is associated with a small energy scale so, equivalently, in the renormalized classical regime the characteristic energy of the fluctuations becomes smaller than the energy associated to temperature ($\hbar\omega<k_{B}T$). Close to $T_{c}$ and to the Widom line, the correlation length $\xi$ grows rapidly (equivalently, there is “critical slowing down” so the characteristic frequency of the system decreases), and hence the thermodynamic fluctuations $\delta N^{2}$ are dominated by the classical contributions ($\chi(q=0,0)\sim\xi^{2}$ in mean-field). By contrast with the quantum contribution to the nonlocal fluctuations $\delta n_{nl}^{2}$ that is negative, the classical (thermal) contribution to the nonlocal fluctuations $\delta n_{nl}^{2}$ is positive. Hence $\delta n_{nl}^{2}$ becomes less negative close to $T_{c}$ or to the Widom line. Returning to the local fluctuations $\delta n^{2}$, the classical (i.e. thermal) contribution to $\delta n^{2}$ is logarithmic in $\xi$, as one can see using the Ornstein-Zernike form for the susceptibility $(\int_{0}^{\Omega}\chi(q,0)d^{2}q\sim\int_{0}^{\Omega}\frac{A}{\xi^{-2}+q^{2}}d^{2}q\sim A\ln(\Omega/\xi))$. This means that $\delta n^{2}$ does not depend much on $\xi$. It is controlled mainly by quantum effects, even though strictly speaking there is a diverging contribution at $T_{c}$.

Therefore there are two mechanisms governing the temperature behavior of $\delta n_{nl}^{2}$ away from half-filling, $n<1$, for a doped Mott insulator. The first mechanism occurs at high temperature, as discussed in Sec. III. As the temperature decreases below the Fermi energy, quantum mechanical effects become important when $\xi_{th}$ becomes larger than the average distance between fermions $\ell$ and it is $\xi_{th}$ then that sets a limit to the range of density fluctuations. The non-local correlations are negative because of the Pauli principle and because of fermion-fermion repulsion. A second mechanism occurs at much lower temperatures: as $T$ crosses the Widom line, or the critical endpoint, the system becomes renormalized classical: $\xi$ becomes larger than $\xi_{th}$ and it is $\xi$ that sets a limit to the range of density fluctuations. Since classical contribution to $\delta n_{nl}^{2}$ are positive, $\delta n_{nl}^{2}$ becomes less negative close to $T_{c}$ or to the Widom line.

V Predictions for experiments

Our results lead to the following predictions for quantum simulations of the two-dimensional Hubbard model with ultracold atoms in optical lattices. For $U$ larger than the critical value $U_{\rm MIT}$ where the Mott transition occurs at half-filling, both the charge compressibility $\kappa_{T}$ and the thermodynamic density fluctuations $\delta N^{2}$ will show a peak as a function of the occupation $n$. As $T$ is decreased toward $T_{c}$, the peak will increase in magnitude and narrow in size. This enhancement in $\delta N^{2}$ is associated with the nonlocal contribution $\delta n_{nl}^{2}$ becoming less negative, while the local fluctuations $\delta n^{2}$ do not do much. These will be the signature of the Widom line. Below $T_{c}$, a first-order transition occurs at finite density and finite temperature, and hence $\delta N^{2}$ will be discontinuous.

Usual provisos apply. First, the presence of long-range ordered phases, or quasi-long-range ordered phases such as antiferromagnetism Mazurenko et al. (2017); Fratino et al. (2017a) or superconductivity Sordi et al. (2012b); Fratino et al. (2016) in the low temperature region of the $T-n$ phase diagram can prevent the observation of the endpoint or of some portion of the Widom line. Nevertheless, the Widom line extends far away from the endpoint, up to $T\approx 5T_{c}$, so it may come out into the normal state. Furthermore, it is possible to shrink the regime of competing states, such antiferromagnetism, by adding frustration (e.g next nearest neighbor hopping $t^{\prime}$ Fratino et al. (2017b), or triangular lattice Hébert et al. (2015)). The onset temperature $T_{N}$ of the antiferromagnetic correlations strongly depends on the interaction $U$ and occupation $n$. Within plaquette CDMFT, it is largest at $n=1$ (i.e. zero doping) around $U\approx 7$, where $T_{N}\approx 0.33$ Fratino et al. (2017a). It decays rapidly as a function of $U$ and $n$. For instance, for $U=6.2$ considered in Sec. IV, antiferromagnetism sets in from $n=1$ to $n\approx 0.85$ Fratino et al. (2017b), similarly to recent experimental data with ultracold atoms in the antiferromagnetic phase Mazurenko et al. (2017). Therefore the results of this work at half filling $n=1$ of Fig. 2 are quantitatively correct down to at least $T\approx 0.33$. Away from half filling, the results for $U=6.2$ in Figs. 3,4 are quantitatively correct for $n\lessapprox 0.85$, but only qualitative for $n\gtrapprox 0.85$.

Second, our results have been obtained with plaquette CDMFT that treats exactly local and short-range fluctuations within the plaquette. The fate of the first order transition and its associated Widom line with larger cluster size including larger range of fluctuations is computationally very challenging and hence has not been explored systematically. Nevertheless, one of our key findings is that nonlocal fluctuations decrease in absolute value near the Widom line and close to $T_{c}$. This suggests that larger nonlocal correlations probed by larger clusters should not qualitatively modify our results.

The significance of the Widom line emanating from the endpoint at finite $T$ and $n$ stems from the theoretical proposal Sordi et al. (2012a) that identifies the Widom line with the rapid crossovers marking the pseudogap temperature $T^{*}$ in the phase diagram of hole-doped cuprates Sordi et al. (2012a, 2013). If confirmed by experiments with ultracold atoms, this result could throw new light into the explanation of the pseudogap phase in cuprate high-temperature superconductors.

Acknowledgements.

This work has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) under grants RGPIN-2014-04584, the Canada First Research Excellence Fund and by the Research Chair in the Theory of Quantum Materials. Simulations were performed on computers provided by the Canadian Foundation for Innovation, the Ministère de l’Éducation des Loisirs et du Sport (Québec), Calcul Québec, and Compute Canada.