  We give a quasi-polynomial time classical algorithm for estimating the ground
state energy and for computing low energy states of quantum impurity models.
Such models describe a bath of free fermions coupled to a small interacting
subsystem called an impurity. The full system consists of $n$ fermionic modes
and has a Hamiltonian $H=H_0+H_{imp}$, where $H_0$ is quadratic in
creation-annihilation operators and $H_{imp}$ is an arbitrary Hamiltonian
acting on a subset of $O(1)$ modes. We show that the ground energy of $H$ can
be approximated with an additive error $2^{-b}$ in time $n^3 \exp{[O(b^3)]}$.
Our algorithm also finds a low energy state that achieves this approximation.
The low energy state is represented as a superposition of $\exp{[O(b^3)]}$
fermionic Gaussian states. To arrive at this result we prove several theorems
concerning exact ground states of impurity models. In particular, we show that
eigenvalues of the ground state covariance matrix decay exponentially with the
exponent depending very mildly on the spectral gap of $H_0$. A key ingredient
of our proof is Zolotarev's rational approximation to the $\sqrt{x}$ function.
We anticipate that our algorithms may be used in hybrid quantum-classical
simulations of strongly correlated materials based on dynamical mean field
theory. We implemented a simplified practical version of our algorithm and
benchmarked it using the single impurity Anderson model.
