  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 $\mathcal{O}(\beta)$ to $\tilde{\mathcal{O}}(\beta^{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\'enyi
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 $e^{\sqrt{\tilde{\mathcal{O}}(\beta \log(n))}}$. This proof 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 $\beta = 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.
