  An expansion method for perturbation of the zero temperature grand canonical
density matrix is introduced. The method achieves quadratically convergent
recursions that yield the response of the zero temperature density matrix upon
variation of the Hamiltonian. The technique allows treatment of embedded
quantum subsystems with a computational cost scaling linearly with the size of
the perturbed region, O(N_pert.), and as O(1) with the total system size. It
also allows direct computation of the density matrix response functions to any
order with linear scaling effort. Energy expressions to 4th order based on only
first and second order density matrix response are given.
