We show that Monte Carlo sampling of the Feynman diagrammatic series (DiagMC) can be used for tackling hard fermionic quantum many-body problems in the thermodynamic limit by presenting accurate results for the repulsive Hubbard model in the correlated Fermi liquid regime. Sampling Feynman's diagrammatic series for the single-particle self-energy we can study moderate values of the on-site repulsion (U/t ∼ 4) and temperatures down to T /t = 1/40. We compare our results with high temperature series expansion and with single-site and cluster dynamical mean-field theory.PACS numbers: 02.70. Ss, 05.10.Ln, 71.10.Fd, 05.30.Fk Advancing first-principle simulations of fermionic manyparticle systems is notoriously difficult due to an exponential computational complexity [1]. One of the most famous examples is the Hubbard model [2, 3] whose phase diagram remains highly controversial,whereĉ † iσ creates a fermion with spin projection σ =↑, ↓ on site i,n iσ =ĉ † iσĉiσ , . . . restricts summation to neighboring lattice sites, and t, U and µ are the hopping amplitude, the on-site repulsion, and the chemical potential respectively.The Hubbard model is of great technological and scientific importance since it is believed by many to be the central model for high-temperature superconductors [3]. Moreover, it can nowadays be realized with ultracold atoms in optical lattices [4] where Mott physics [5,6] has been observed recently. The prospect of using cold atomic gases as building blocks of quantum simulators is thus becoming a reality. However, claiming the equation of state, either numerically or experimentally, can only be substantiated on the basis of unbiased tools. Not only numerical methods need to be tested against known unbiased results, but also quantum emulators need to be validated by benchmarking them at ever lower temperatures. In the bosonic case this can be done with high precision [7], but for fermions it is still an open question what method can be used to accurately describe the lowtemperature regime of the model. Precise numerical solutions of the Hubbard model can only be obtained in the particlehole symmetric case n σ = 1/2, where the sign problem is absent in quantum Monte Carlo (QMC) simulations. In all other cases circumventing exponential scaling necessarily involves approximations with systematic errors that are hard to assess and control.In general terms, Diagrammatic Monte Carlo (DiagMC) [8] is a numerical approach based on stochastic sampling of Feynman diagrams [9, 10] to increasingly higher orders in the coupling constant-rather than evaluating integrals over internal variables for each diagram, one samples their momenta, imaginary times and expansion orders stochastically. So far, DiagMC algorithms have been successfully used for obtaining unbiased solutions for polaron problems (for the latest work, see Ref.[11]).Nevertheless, the crucial question remained unclear of whether the DiagMC approach could be applied to true manybody systems in the thermodynamic limit at physically interesting temper...