{'fraction_non_alphanumeric': 0.04152184871720039, 'fraction_numerical': 0.015627072775994268, 'mean_word_length': 4.175339832486612, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
N\'eel temperature and thermodynamics of the half-filled 3D Hubbard model by Diagrammatic Determinant Monte Carlo

We study thermodynamics of the 3D Hubbard model at half filling on approach to the N\'eel transition by means of large-scale unbiased Diagrammatic Determinant Monte Carlo simulations. We obtain the transition temperature in the strongly correlated regime, as well as temperature dependence of energy, entropy, double occupancy, and the nearest-neighbor spin correlation function. Our results improve the accuracy of previous unbiased studies and present accurate benchmarks in the ongoing effort to realize the antiferromagnetic state of matter with ultracold atoms in optical lattices.


I. INTRODUCTION
The Hubbard model 1 of interacting fermions in a solid is a centerpiece of modern condensed matter physics. It is conventionally defined by restricting the motion of electrons in a crystalline solid to a single band, and simplifying the screened long-range Coulomb interactions between electrons to an on-site repulsion: (1) where σ =↑, ↓, c † xσ creates a fermion on a site x, n xσ = c † xσ c xσ , the summation in the first term runs over the nearest-neighbor sites of the simple cubic lattice, t is the hopping amplitude and U > 0 is the onsite repulsion. Though remarkably simple in appearance, the model has been used to study a wealth of intriguing quantum many-body phenomena that are due to electron correlations in solids, such as interaction-driven insulators 2 , quantum magnetism and high-temperature superconductivity 3 . However, despite more than a half of century of intensive investigation, the physics of the model is still not completely understood.
The most challenging yet the most interesting regime is the intermediate regime with interaction comparable to half the bandwidth U ∼ zt, where z is the number of nearest neighbors of a site. This regime offers no small parameter to start a controllable analytic theory. Furthermore, exact analytic solutions are only accessible either in one spatial dimension 4 or an infinite 5 number of spatial dimensions. Substantial progress has been possible with the development of efficient quantum Monte Carlo (QMC) methods (for a recent review, see Ref. 6) accompanied by advances in computer technology. Although for generic bosonic systems virtually any equilibrium property can nowadays be calculated by QMC with a controlled high accuracy 7 , systematic-error-free simulations of correlated fermions have been limited to a handful of special cases due to the negative sign problem 8 . The sign problem manifests itself as an exponential scal-ing of the simulation time with the system volume and inverse temperature making it practically impossible to obtain any reliable information about the system in the thermodynamic (TD) limit. Although in some cases the sign problem can be completely eliminated by choosing a system-specific representation, its general solution is almost certainly not possible 8 .
A major step toward understanding strongly correlated systems has been the experimental realization of the Hubbard model with ultracold atomic gases loaded into optical lattices (for recent reviews see Refs. 9 and 10). These systems offer substantial control over the Hamiltonian. As a result, these experiments can serve as emulators of quantum many-body systems, which allow the accurate study of a given model in a range of parameters inaccessible by analytic and numeric techniques 11 . The recent experimental observation of Mott physics 12,13 in the Hubbard model is a major milestone along these lines. The next crucial step would be a realization of the antiferromagnetic (AFM) transition and the Néel state in the Hubbard model, which requires a substantial effort in reaching lower temperatures, controlling the equilibration rates, developing new probing techniques, etc. In addition to these inherent challenges, there is also a fundamental problem related to thermometry in ultracoldatom systems since they are insulated from the environment. Isolated systems such as these are thus, by their nature, better characterized by entropy, rather than temperature. Moreover, probes have to be calibrated in the relevant regime and the results obtained with the setup validated against available benchmarks. For these purposes, reliable and accurate numeric results for fermionic systems are indispensable because they ultimately allow a full quantitative understanding, as was recently demonstrated by the example of a bosonic optical-lattice emulator 7 .
This work provides reliable benchmarks for the realization of the Néel state in optical lattices as well as for new theoretical methods. We focus on the special case of the half-filled ( n xσ = 1/2, or equivalently µ = U/2) three dimensional (3D) Hubbard model (1) on a simple arXiv:1212.3027v2 [cond-mat.quant-gas] 24 May 2013 cubic lattice. The case of half filling is special due to the SU (2) symmetry of the Hamiltonian, which is ultimately broken by the Néel state, making magnetism the leading instability. In the limit U/t 1, the effects of the interaction can be studied perturbatively. In the opposite limit, U/t 1, Eq. (1) reduces to the AFM Heisenberg model with J ∼ t 2 /U . While there is no doubt that in the strongly correlated regime the ground state of the half-filled model (1) is AFM, mapping out the finitetemperature phase diagram and studying thermodynamics of the system is extremely challenging [14][15][16][17][18] . To this end, we employ the unbiased continuous-time determinantal diagrammatic Monte Carlo (DDMC) technique 19 , which produces numerically exact (up to a known statistical error bar) results for a finite-size system, and which is free from the fermionic sign problem at half filling on bipartite lattices, allowing a reliable extrapolation of results to the thermodynamic limit.
We use DDMC to determine several critical properties of the Hubbard model with high control and accuracy and compare with high-temperature series expansions 20 (HTSE) where possible. We study the range of on-site interaction 4 ≤ U/t ≤ 8, where the critical temperature T N of the AFM transition is expected to reach its maximum 14 . Results in this regime are vital to optical lattice emulator efforts offering experiments their best chance of observing the AFM phase. We obtain the critical temperature T N and compute important thermodynamic properties of the model-the energy and the entropy, as well as two optical lattice observables: double occupancy and the nearest-neighbor spin-spin correlation function-in the paramagnetic phase as a function of temperature down to T N . A number of previous works using different unbiased approaches studied T N (by the determinant quantum Monte Carlo (DQMC) 14 and the Dynamical Cluster Approximation (DCA) 15 ) and the thermodynamic properties in question (by DDMC 16 , a combination of DCA and DDMC 17 , and by DQMC 18 ) in this regime. Our work improves the accuracy of the previous results at half filling and provides the most accurate to date values of the critical temperature T N and the entropy at the critical point S N (summarized in the Table. I). Accurate knowledge of T N is particularly important for determining the critical entropy S N since S(T ) is a steep function near the transition, so that the error bar of S N is mainly due to the uncertainty in T N . The values of S N are required for an experimental realization for the AFM state. As was noted in Ref. 17, close to the transition, the temperature dependence of the nearestneighbor spin-spin correlation function S z x S z x+ei (e i is the unit vector in the direction i) is significantly more pronounced than that of the double occupancy n x↑ n x↓ of a lattice site making measurements of the nearestneighbor spin-spin correlations 21,22 more suitable for accurate thermometry in this regime. Our results for the spin-spin correlation function can be used for calibration of such a thermometer.
For over a decade, the DQMC study of the phase di-agram of the half-filled 3D Hubbard model by Staudt et al. 14 has been the main reference for T N in the correlated regime. Representing the state of the art at that time, Ref. 14 provides a comprehensive comparison of the DQMC data for the Néel temperature with those of preceding QMC simulations and approximate theories, e.g., DMFT 23 . We shall not reproduce this comparison here and refer the reader to a review of results for T N predating Ref. 14.
The simulation method of Staudt et al. is based on a discrete Hubbard-Stratonovich decoupling in the Hubbard interaction term which requires discretizing the imaginary-time interval 0 < τ < β = 1/T into a finite number of steps of size ∆τ thereby introducing a systematic error. Hence, in addition to the standard extrapolation of the finite-size DQMC data to the TD limit, one has to perform an extrapolation with respect to ∆τ → 0. Such a double extrapolation is rather laborious. In practice ∆τ is usually fixed at a value which is large enough to allow efficient simulations, yet, according to Ref. 14, such that "the results are not significantly affected by the extrapolation ∆τ → 0". In the absence of an explicit extrapolation, the degree of control over systematic errors can be questioned. This is where our approach is a major improvement over that of Staudt et  We have also been able to identify the transition itself with a considerable improvement in accuracy over Ref. 14. Long-range AFM order in the system causes divergence of the magnetic structure factor where s z (x) = (n x↑ − n x↓ )/2, and Q = (π, π, π) is the AFM wave vector, so that S(Q)/L 3 is related to the magnetization m in the TD limit, lim L→∞ S(Q)/L 3 = m 2 . In Ref. 14, the transition temperature is found as the point at which m starts to noticeably deviate from zero, while the magnetization is obtained from a finite-size extrapolation of S(Q) with respect to L → ∞. Here the scaling power of the finite-size corrections was used as a fitting parameter. An additional (indirect) probe of the transition used by Staudt et al. is the peak in the dependence of the specific heat on temperature. In contrast, we use a much more accurate method to determine the critical point. We find the critical point using a finite-size scaling analysis of S(Q) in combination with the technique of Binder crossings 33 . This approach allows us to get a reliable and accurate value of T N by making use of the known scaling law of the magnetic structure factor at the critical point, S(Q) ∝ 1/L −2+η , where η is the anomalous dimension of the particular universality class, so that the quantity S(Q)L −2+η becomes scale-invariant at the critical temperature up to non-universal corrections, which we also take into account. Since the Néel transition is breaking the SU (2) symmetry of the Hubbard Hamiltonian at half-filling, the universality class is that of the 3D Heisenberg model, which provides the critical exponents for the finite-size scaling analysis.
In the context of previous calculations of the Néel temperature, it is worth noting the work by Kent et al. 15 , where T N was found from calculations using DCA. Although this work did not lead to any improvement of the precision claimed by Staudt et al., it was shown that the DCA can be used to determine the critical point with a controlled accuracy based on significantly smaller clusters of only ∼ 50 lattice sites versus up to ∼ 1000 sites in Ref. 14 and in this work. In addition, good agreement with the values of T N of Ref. 14 suggested that potential systematic error of the ∆τ → 0 extrapolation in Ref. 14 is likely to be small. In the range of interaction U/t ≥ 6, our results for T N agree within the error bars with those of Refs. 14 and 15 implicitly confirming this. At smaller U , however, we find somewhat lower values of T N . Moreover, at U/t = 4, being unable to reach significantly low temperatures to accurately infer the critical point, we can only claim an upper bound from our finitesize-scaling analysis, T N /t < 0.17, which is already lower than the values claimed in Refs. 14 and 15. The reason for the discrepancy is likely to be the long-range character of correlations at smaller coupling, which can be missed in simple finite-size extrapolation schemes based on data for insufficiently large systems.
The thermodynamics of the Hubbard model near the Néel transition in connection with its experimental realization has been the focus of a number of recent studies [16][17][18][24][25][26][27][28][29] . The DMFT results 26,29 emphasize the role of double occupancy in detecting the build up of AFM correlations. However, its dependence on temperature near the Néel transition is relatively flat in the regime where T N is maximal, as observed in unbiased (extrapolated to the TD limit) DCA calculations 17 . Fuchs et al. 17 obtained the energy and the entropy down to T N in this regime as well as the equation of state away from half filling, which allowed to get an estimate of the entropy at the transition in a realistic harmonically trapped system. A study 18 of the system using the same DQMC method as Ref. 14 arrived at an agreement with the results and conclusions of Fuchs et al..
The DDMC simulation method used in this work is not capable of capturing thermodynamics of the Hubbard model away from half filling due to the pronounced negative sign problem. However, exactly at half filling, it has certain advantages over DCA and DQMC. As discussed above, in order to claim unbiased results in the TD limit within DQMC one has to resort to a double extrapolation, ∆τ → 0, L → ∞. DDMC is formulated directly in continuous time allowing a more reliable extrapolation to the TD limit. Modern efficient solvers for DCA also work in continuous time, but the cluster sizes amenable to simulation in the regime of interest are typically less then 100. However, in DCA the clusters are embedded in a self-consistently defined medium, which largely improves the convergence to the TD limit. In practice the finite-size dependence for the accessible clusters is notably larger than that in DDMC 17 . Hence, at half filling and U/t ∼ 8, where T N and S N are expected to reach their maxima, DDMC currently allows to obtain the most reliable benchmarks for the thermodynamic quantities of interest.
The paper is organized as follows. In Sec. II, we discuss the simulation method outlining the general formulation of the DDMC technique and its application to calculating the specific observables in question. Section III is concerned with determining the temperature of the Néel transition. Sec. IV describes the thermodynamics near T N . Here we discuss the extrapolation of the observables to the TD limit (IV A), the determination of entropy (IV C), and thermometry near T N (IV D). We summarize the results in Sec. V. The Appendix contains tables of the obtained numerical data for the entropy, energy, double occupancy, and spin-spin correlation functions as a function of temperature.

II. METHOD
We first rewrite the Hubbard Hamiltonian (1) to a form suitable for numerical simulations by mapping the repulsive model (1) to an attractive model by a particle-hole transformation 30 . We use the fact that the simple cubic lattice is bipartite and can be split into two interpenetrating sublattices A and B, so that the hopping term in (1) only connects sites belonging to different sublattices. Then we introduce the hole operators for the ↑component: This way, Eq. (1) becomes at half filling where a x↓ = c x↓ , m xσ = a † xσ a xσ is the number operator for the attractive model, µ = −U/2 as appropriate for half filling, and Ω = L 3 is the total number of sites.
Since we only consider half filling, m x↑ = m x↓ = 1/2, we follow Ref. 31 and shift the chemical potential according to µ → µ + αU : with At half filling the choice α = 1/2 is optimal because it leads to the minimal computational complexity of the simulations (see below), and we use this value of α throughout.

A. Diagrammatic Determinantal Monte Carlo method
The DDMC algorithm works with the weak-coupling expansion series for the finite-temperature partition function for the Hubbard model (6). The latter reads where Here β is the inverse temperature, Z 0 = Tr T exp (−βH 0 ) is the unperturbed partition function, T denotes time ordering, and (. . . ) 0 = Tr [T (. . . ) exp (−βH 0 )] /Z 0 denotes the thermodynamic average with respect to the unperturbed Hamiltonian H 0 . Equations (7)-(8) generate the standard Feynman diagrams: There are (p!) 2 diagrams of order p, which can be represented graphically as a collection of p vertices connected by single-particle propagators. Summing over all the interconnections for a fixed vertex configuration equation (8) takes the form where A(S p ) are p × p matrices with matrix elements given by (i, j = 1, . . . , p) (since we only consider an unpolarized system, we omit the spin index for A-s and G (0) -s), and G (0) being the free-particle Green's functions, SinceĤ 0 , Eq. (6), is diagonal in momentum space, the free propagators are calculated by the Fourier transform , τ > 0, (14) and tabulated before the start of the simulation. The series (7), (10)-(11) serves as a basis for a DDMC simulation of the Hubbard model (6): We set up a random walk in the space of the vertex configurations S p , Eq. (9), using the Metropolis algorithm 32 with the weights proportional to D(S p ), Eq. (10). Since the technique itself is detailed elsewhere 19,31 , here we only briefly discuss the specific details of the present implementation. We only stress at this point that since all the terms in the series (7),(10)-(11) are positive definite we completely avoid a sign problem.
The simplest updating strategy for DDMC simulations consists of adding (p → p + 1) and removing (p → p − 1) interaction vertices at random positions in x and τ to/from a vertex configuration S p . However, at half filling and with α = 1/2 the series only contains even order terms, hence we employ rank-2 updates, where p → p ± 2. Using the Woodbury-type formulas, both rank-1 and rank-2 updates can be performed in O(p 2 ) operations. We note that for α = 1/2 both even and odd terms are present even at half filling. In this sense the choice of α = 1/2 is optimal.

B. Observables
The general method for calculating observables in the DDMC simulations uses the standard technique of Monte Carlo estimators: for an observable O we define an estimator Q (O) (S p ) such that the average of the latter over the vertex configurations generated by the MC process, Below we explicitly list estimators for useful observables: a. Filling fraction The thermal average of the filling fraction m σ for the spin projection σ is given by hence the corresponding estimator is 19 where A(S p ) is a p × p matrix (11) and B(S p ; xτ ) is a (p + 1) × (p + 1) matrix with an extra row and a column corresponding to the extra creation and annihilation operators in (16). Here x and τ are random positions in space and time, respectively. Notice that for a half filled model and for α = 1/2, the MC average of the second term in (17) must equal zero. We use this fact to check if a simulation has equilibrated.
b. Kinetic energy Calculating the kinetic energy for the model (6) requires evaluating the average of a † xσ a yσ , which differs from Eq. (16) only in that the creation annihilation is shifted in space with respect to the creation operator while in Eq. (16) both operators reside on the same lattice site. The estimator for the kinetic energy H 0 is then 19 where the matrix B(S p ; x, y, τ ) differs from B(S p ; xτ ) by the fact that the creation operator in (18) is shifted in space with respect to the annihilation operator. The extra factor of 2 accounts for the summation over σ =↑, ↓ and dL 3 is the number of bonds of a lattice with L 3 sites and periodic boundary conditions (PBC). c. Interaction energy Since the series (7) is nothing but an expansion in powers of H 1 , the corresponding estimator is readily obtained by a standard trick of considering the Hamiltonian H 0 + λH 1 and differentiating with respect to λ . The result is 19 d. AFM structure factor Calculating the AFM structure factor (2) in the particle-hole transformed model (6) requires calculating two independent equaltime density-density correlation functions The estimator for the equal spin density correlation function g ↑↑ for x = y is given by where B 2 (S p ; x, y, τ ) is a (p + 2) × (p + 2) matrix with two extra rows and columns corresponding to the extra creation and annihilation operators in Eq. (20). For x = y, g ↑↑ (x = y) = m 2 x↑ which equals 1/2 for a half-filled model.
To build the estimator for Eq. (21) we proceed similarly to (19). The resulting estimator for (21) is where B(S p ; y, τ ) is a p × p matrix constructed by selecting a random vertex, (xτ ), from a configuration S p and moving the corresponding row of the matrix A(S p ) to (yτ ) while leaving the corresponding column at (xτ ) .

III. CRITICAL TEMPERATURE
In the paramagnetic phase, T > T N , S(Q) scales to zero exponentially as L → ∞. In the AFM phase, on the other hand, S(Q)/L 3 → m 2 as L → ∞, where m is the sublattice magnetization. Right at the critical temperature S(Q) ∝ 1/L −2+η , where η is the anomalous dimension.
In order to locate the transition temperature we thus use the standard finite-size scaling (FSS) ansatz 33 where ξ is the correlation length which diverges at the transition as ξ ∝ |T − T N | −ν , f (x) is a real-valued function which tends to a finite constant as x → 0, and the corrections in brackets arise from the leading irrelevant operators (dots represent the higher-order corrections).
Here the exponent ω is universal, but the amplitude c is not. In accordance with the 3D Heisenberg universality class, we take η ≈ 0.037, ν ≈ 0.71 and ω ≈ 0.8 34 .
The basic idea of using Eq. (24) for the FSS is as follows: if the corrections-to-scaling (the 2nd term in brackets in Eq. (24)) were not present, S(Q)L −2+η would be scale independent at the transition point, so that performing the simulations at a series of system sizes L 1 > L 2 > . . . and plotting S(Q)L −2+η versus temperature, one would observe that all the curves intersect at the same point, T = T N . This is what we indeed observe for (relatively) large values of U/t at L ≥ 6: for U ≥ 6t our MC results are consistent with c = 0 in Eq. (24) within statistical errors, see Fig. 1. The data for the smallest systems of size L = 4 systematically deviate from the scaling described by Eq. (24) for all values of U considered here (see also below), and hence they are omitted from the scaling analysis.
We find that the corrections-to-scaling become more pronounced with decreasing U/t: For U = 5t we find a clear evidence of the shift of the pairwise crossings towards lower temperatures, see Fig. 2. Since simulating larger system sizes is not an option, we employ Eq. (24) including corrections to scaling. The most straightforward way is to follow the evolution of the pairwise crossings with the system size. Indeed, expanding Eq. (24) around the crossing of S(Q)L −2+η at system sizes L = L 1 and L = L 2 up to the terms linear in T − T N we find (cf Ref. 19):  where g(L 1 , L 2 ) = 1 We perform the linear fit of the series of crossings T L1,L2 versus g(L 1 , L 2 ). Then the intercept of the best-fit line yields the Néel temperature. This procedure is illustrated in Fig. 3. It is clear from Figs. 1-3 that L = 4 does not follow the scaling (25)- (26). We attribute it to the effect of the higher-order terms neglected in (24), and only use L > 4 in the fitting procedure (25)- (26). We note at this point that the methodology based on Binder crossings has a build-in self-consistency check: if, in fact, the criticality were not in the Heisenberg universality class, the curves for the magnetic structure factor, rescaled via Eq. (24) would have no reason to cross at a unique point, and thus the whole procedure of (24)-(26) would break down.
An equivalent procedure has been suggested in Ref. 35. Again, one expands Eq. (24) up to the linear order in T − T N , which leads to which is then used as a four-parameter ansatz for a single nonlinear fit. A priori, fitting procedures based on (25) and (27) are equivalent and indeed produce consistent results, as illustrated in Fig. 2. We stress at this point that using (27) requires judicious choice of the temperature range for fitting: including Monte Carlo points at too high temperatures and/or too small system sizes tends to significantly skew the fit results. In the following we therefore quote the T N -s obtained using Eq. (25).
For U = 4t, we find the corrections-to-scaling to be larger than those for U = 5t, see accessible systems sizes we are only able to put an upper limit on the Néel temperature, T N < 0.17t. From Fig. 4 it is clear that L = 6 and possibly even L = 8 are simply too small and need to be discarded from the finite-size scaling analysis.
Our results for the dependence of the Néel tempera-   ture on U are summarized in Table I and Fig. 5. It is instructive to compare our estimates to the previous unbiased calculations from the literature. While for U/t = 6 and 8 our estimates agree with and are more accurate than previous estimates from QMC 14 and DCA 15 . For smaller values of U/t our estimates are systematically lower. The discrepancy can be traced back to the FSS procedure which includes corrections-to-scaling, Eq. (24): if we were to discard the corrections and identified the Binder crossings of L = 6 and L = 8 with the Néel temperature, such estimates would have agreed with Refs. 14 and 15. We therefore conclude that the estimates of T N presented here are more accurate than results reported to date.

A. Extrapolation to the thermodynamic limit
The dependence of local observables on the size L of the system with PBC is complicated 36 by oscillations between the results for even and odd values of L/2. The issue is illustrated in Fig. 6, where the energy per particle  of the half-filled non-interacting system (U = 0, µ = 0) with PBC is plotted versus L −1 up to a large system size (L = 52) for different temperatures. The TD-limit value is approached from above by the data for even L/2 and from below by those with L/2 odd. These are the well-known "shell" oscillations 37 caused by whether or not the spectrum of the finite system has states ν with the energy E ν within a range much less than T from the Fermi level, |E ν − µ| T . In the example of Fig. 6, the states are classified by the momenta k = (k 1 , k 2 , k 3 ), k i = 2πn i /L with integers n i taking the values n i = −L/2, . . . , −1, 0, 1, . . . , L/2 − 1. When {n i } = L/4, which is only possible if L/2 is even, the state k = (π/2, π/2, π/2) is exactly at the Fermi level; it's occupation is 1/2 ("open shell"), but it gives no contribution to the average energy. Hence, if L is not large enough so that the spacing between the levels is larger than T , the average energy per particle of the system with a closed shell (L/2 odd) is systematically lower than that of the system with an open shell (L/2 even) due to the difference of the number of states below the Fermi level, L 3 /2 and L 3 /2 − 1 correspondingly. However, for any given temperature T there is a system size L * = L * (T ) such that for L > L * the number of states with the energies |E ν − µ| < T becomes large removing the distinction between even and odd L/2. In the free-particle case of Fig. 6, the convergence to the TD limit at L > L * is extremely fast (exponential) with L * (T = 0.5) ∼ 10, L * (T = 0.3) ∼ 16, and L * (T = 0.1) ∼ 40.
At the finite values of U studied here, the paramagnetic phase should be described by a Fermi liquid in the limit of T E F . In this regime, the total energy is a functional of occupation numbers of noninteracting quasiparticles. Therefore, the system-size dependence of energy is expected to be proportional to that of the non-interacting system, at least for large enough   Fig. 7. The figure suggests that the additional corrections given by g(L) should be small for L ≥ 6 at large U , whereas they are appreciable for most of the accessible system sizes at smaller U . We claim the TD-limit value E(L → ∞) and its error bar ∆E(L → ∞) conservatively as the span between the values at the two largest accessible system sizes including their statistical error bars (depicted by the horizontal band in the upper panel of Fig. 7): and and similarly for other local observables. As a consistency check for our TD-limit results (at U ≥ 6) as well as to improve convergence to the TD limit (at U ≤ 5), we employ two other simulation setups, which exhibit different system-size dependences, which we detail below.  (9) for U = 4, T = 0.2 (using the TABC data)-are shown by the horizontal bands. More generally, we use PBC data at U ≥ 6 and TABC data at U ≤ 5 to obtain the thermodynamic-limit values, as explained in Subsection. IV A.

Twist-averaged boundary conditions.
Averaging over twisted boundary conditions was found in Refs. 37 and 39 to produce exact results for the noninteracting system in the grand canonical ensemble and to substantially suppress the system-size dependence for interacting systems. In this approach one introduces a finite phase that particles acquire when they wrap around the periodic boundaries, |r 1 + Le i , r 2 , . . . = e iΘi |r 1 , r 2 , . . . , i = 1, 2, 3, (30) where e i is the unit vector in the direction i and −π < Θ i ≤ π; Eq. (30) with Θ i = 0 corresponds to the standard PBC. Then an observable A(L) is obtained by means of the integration where A Θ (L) is a result of the simulation with a fixed value of Θ = {Θ i } and system size L. Thereby, the possible momentum values are forced to span the whole Brillouin zone. The non-interacting propagators G (0) satisfying the condition (30) are obtained by substituting k → k + Θ/L in Eq. (14). In practice, we perform numerical integration on a mesh of 64 Θ points, estimating the systematic error of integration to be smaller than the statistical error of A(L) coming from sampling each A Θ (L) by Monte Carlo.
The results of the calculation of energy with the twistaveraged boundary conditions (TABC) are compared to those for the PBC in Fig. 7.
For U ≤ 5 we find the TABC to substantially reduce the finite size corrections, as exemplified by the lower panel of Fig. 7 showing the typical comparison data for U = 4. Correspondingly, at U ≤ 5 we base our L → ∞ extrapolation on the TABC data and use Eqs. (28), (29) with E[L] and ∆E[L] being the finite-size value and its error bar obtained with TABC to claim the TD-limit extrapolated results (exemplified by the horizontal band in the lower panel of Fig. 7).
For larger values of U/t, averaging over the twists still reduces the finite size corrections somewhat, but the net improvement of TABC over PBC is smaller (for a typical example see the upper panel of 7-notice that for U = 8 data for L = 8 and L = 10 appear to be converged within their error bars to a value consistent with the PBC extrapolated value) Overall, we find that for U/t ≥ 6 the use of TABC does not lead to a significant improvement of the convergence to the TD limit, and we thus use PBC data and Eqs. (28), (29) in this range of interactions. We follow the same protocol to obtain other observables.

L → ∞ free propagators.
In the second approach, we replace the free-particle propagators G (0) = G ∞ , thereby completely eliminating the oscillations coming from the discreteness of the spectrum at the expense of giving up the PBC. The only source of systematic error in this case is the finiteness of the volume-still given by L 3 -confining the distribution of the interaction vertices in Eq. (9). In this case, the finite-size corrections are substantially larger than those of simulations with PBC. However, the scaling of these corrections is linear in 1/L for all the local observables in question, which allows to perform a systematic TD-limit extrapolation. As an example, such an extrapolation for energy in comparison with the data for PBC (at U = 8) and TABC (at U = 4) is shown in Fig. 8. The TD-limit value obtained thereby is in perfect agreement with the result of simulations with the PBC and TABC. This constitutes an independent verification of the accuracy of the claimed results.

B. Energy
Here we present simulation results for the total energy per particle extrapolated to the TD limit for a range of the interaction U in the correlated regime near the Néel transition. The temperature dependence of the energy per particle for U = 8, 6, 5, 4 is plotted in Fig. 9 along with the results of the HTSE 20 for orders 2, 4, 6, 8, 10. As seen from the plot, the HTSE starts diverging well above the transition point.

C. Entropy
The entropy per particle S(T ) at a given temperature T is obtained from the thermodynamic relation T dS = dE at fixed volume by the integral where T * is some temperature at which the entropy is known. We choose T * to be the lowest temperature at which the HTSE for the energy obviously converges to the TD-limit value E(T * ) from the simulation. From Fig. 9, we find T * = 1.8, 2.4, 2.6, 2.6 at U = 8, 6, 5, 4 respectively. Then, the accurate value of S(T * ) in Eq. (32) is given by the HTSE, while the integral is done over the simulation data after taking the TD limit. Since the dependence E(T ) is slow, we represent it by a piecewise linear function and take the integral analytically. The systematic error of integration is included in the error bars for S(T ), but is negligible compared to the error propagated from the values of E(T ). The resulting curves of S(T ) for U = 8, 6, 5, 4 are shown in Fig. 10. From these data and our calculation of T N discussed in Sec. III, we find the values of the critical entropy S N = S(T N ) in the range of U and summarize the results in Table I. The error bars of S N are dominated by the relatively small error of T N due to the large slope of S N (T ) near the transition. In Fig. 11, we plot lines of constant entropy in the (T , U ) plane. The latter demonstrate that an adiabatic increase of the coupling U can lead to either a rise (at S 0. 35  ∞ deviates from the linear scaling followed by larger systems and therefore is excluded from the fit. a fall (at 0.35 S 0.7) of temperature, although the net effect of the Pomeranchuk cooling near S N is rather small.

D. Thermometry
As was shown in Ref. 17 by means of DCA calculations, the nearest-neighbor spin-spin correlation func-tion defined as S z x S z x+ei , which is accessible in presentday ultracold-atom experiments, can serve as a sensitive thermometer at temperatures near T N . In contrast, another routinely measured correlator, the double occupancy n x↑ n x↓ of a lattice site, is nearly flat in this temperature range making it a rather poor candidate for thermometry. This is hardly surprising since the latter is concerned with correlations in the charge channel, whereas the relevant physics at these temperatures is that of developing short-range spin correlations.
In Fig. 12, we present our results for S z x S z x+ei and n x↑ n x↓ extrapolated to the TD limit at U = 8, 6, 5, 4. The obtained values agree within the errors with the TD-limit-extrapolated data from the DCA simulations, Ref. 17, but our error bars are notably smaller. Our data can be directly used for thermometry calibration and detection of the Néel transition. Note that in the range of temperatures T N < T < 2T N (the position of T N is depicted by a vertical line with the width corresponding to the error bar) the spin-spin correlations between nearest-neighbor cites rise by a factor of two, with a substantial increase of the slope close to T N . In the same temperature range the double occupancy varies by less than 10%.

V. CONCLUSIONS
We presented unbiased results for the 3D Hubbard model at half filling near the Néel transition in the strongly correlated regime of 4 ≤ U ≤ 8, where T N reaches its maximum. We focused on the properties of the model near the transition accurately determining T N and studying the energy, entropy, double occupancy n x↑ n x↓ and the nearest-neighbor spin-spin correlator S z x S z x+ei as functions of temperature and interaction. Accurate quantitative understanding of the model in this regime is of growing importance in view of the ongoing experimental effort to emulate the Hubbard model with ultracold atoms in optical lattices, which could ultimately allow to study regions of the phase diagram inaccessible by unbiased theoretical methods. In particular, this could lead to answers of fundamental questions regarding the nature of superfluidity at finite doping and its connection to high-temperature superconductors 11 . The realization of the Néel state would be a necessary step on the way to accessing the region of the phase diagram where quantum fluctuations play an important role. Our simulations provide the most accurate and controlled estimates of entropy at the critical point to date. These entropies, summarized in Table I, have to be achieved in the middle of the trapped cold-atom system to realize the AFM state. For independent in situ thermometry in this regime, one can employ measurements of the nearest-neighbor spin correlations, which expectedly have pronounced temperature dependence near T N , and which can nowadays be addressed either by the use of superlattices 21 or by lattice modulation 22 (3)) is below that claimed in Ref. 17 (0.42(2)) suggesting agreement at the level of two combined standard deviations.