  In one-dimensional quantum systems with short-range interactions, a set of
leading numerical methods is based on matrix product states, whose bond
dimension determines the amount of computational resources required by these
methods. We prove that a thermal state at constant inverse temperature $\beta$
has a matrix product representation with bond dimension $e^{\tilde
O(\sqrt{\beta\log(1/\epsilon)})}$ such that all local properties are
approximated to accuracy $\epsilon$. This justifies the common practice of
using a constant bond dimension in the numerical simulation of thermal
properties.
