One of the most fundamental problems in quantum many-body physics is the characterization of correlations among thermal states. Of particular relevance is the thermal area law, which justifies the tensor network approximations to thermal states with a bond dimension growing polynomially with the system size. In the regime of sufficiently low temperatures, which is particularly important for practical applications, the existing techniques do not yield optimal bounds. Here, we propose a new thermal area law that holds for generic many-body systems on lattices. We improve the temperature dependence from the original O(β) to O~(β2/3), thereby suggesting diffusive propagation of entanglement by imaginary time evolution. This qualitatively differs from the real-time evolution which usually induces linear growth of entanglement. We also prove analogous bounds for the Rényi entanglement of purification and the entanglement of formation. Our analysis is based on a polynomial approximation to the exponential function which provides a relationship between the imaginary-time evolution and random walks. Moreover, for one-dimensional (1D) systems with n spins, we prove that the Gibbs state is well-approximated by a matrix product operator with a sublinear bond dimension of eO~(βlog(n))√. This allows us to rigorously establish, for the first time, a quasi-linear time classical algorithm for constructing an MPS representation of 1D quantum Gibbs states at arbitrary temperatures of β=o(log(n)). Our new technical ingredient is a block decomposition of the Gibbs state, that bears resemblance to the decomposition of real-time evolution given by Haah et al., FOCS'18.
