We show how the thermodynamic properties of large many-body localized systems can be studied using quantum Monte Carlo simulations. To this end we devise a heuristic way of constructing local integrals of motion of very high quality, which are added to the Hamiltonian in conjunction with Lagrange multipliers. The ground state simulation of the shifted Hamiltonian corresponds to a high-energy state of the original Hamiltonian in case of exactly known local integrals of motion. We can show that the inevitable mixing between eigenstates as a consequence of non-perfect integrals of motion is weak enough such that the characteristics of many-body localized systems are not averaged out in our approach, unlike the standard ensembles of statistical mechanics. Our method paves the way to study higher dimensions and indicates that a full many-body localized phase in 2d, where (nearly) all eigenstates are localized, is likely to exist.Introduction -Many-body localization (MBL) addresses the fundamental question under which conditions quantum systems can avoid ergodicity and thermalization, thereby generalizing Anderson localization to interacting systems [1][2][3][4][5][6][7].The widely accepted mechanism for thermalization in quantum systems is the eigenstate thermalization hypothesis [8][9][10]: Under very mild assumptions, (almost) every eigenstate of the system is thermal. This implies that the reduced density matrix of a small subsystem, obtained by tracing out the degrees of freedom of the considered (eigen)state outside the subsystem, is indistinguishable from the thermal density matrix with an effective temperature that depends on the energy density of the chosen eigenstate. MBL states, on the contrary, retain knowledge of their initial local conditions in local operators for asymptotically large times. The picture of local integrals of motion (LIOM) [11][12][13][14][15][16] can explain most of the unusual phenomenology of MBL states: an area low holds in space leading to a logarithmic growth of the entanglement entropy in time and space, and the dc conductivity is identically zero. The area law was demonstrated explicitly in Ref. [17], while it was also found that rare regions can lead to deviations. Dynamics is a decisive characteristic to distinguish between thermal and MBL states: experiments on trapped ions [18] and 1d cold atoms [19] demonstrated memory of the initial conditions over long periods of time for sufficiently strong disorder.