  While general quantum many-body systems require exponential resources to be
simulated on a classical computer, systems of non-interacting fermions can be
simulated exactly using polynomially scaling resources. Such systems may be of
interest in their own right, but also occur as effective models in numerical
methods for interacting systems, such as Hartree-Fock, density functional
theory, and many others. Often it is desirable to solve systems of many
thousand constituent particles, rendering these simulations computationally
costly despite their polynomial scaling. We demonstrate how this scaling can be
improved by adapting methods based on matrix product states, which have been
enormously successful for low-dimensional interacting quantum systems, to the
case of free fermions. Compared to the case of interacting systems, our methods
achieve an exponential speedup in the entanglement entropy of the state. We
demonstrate their use to solve systems of up to one million sites with an
effective MPS bond dimension of 10^15.
