The recursive polynomial expansion for construction of a density matrix approximation with rigorous error control [J. Chem. Phys... 128, 074106 (2008)] is implemented in the quantum chemistry program Ergo [SoftwareX 7, 107 (2018)] using the Chunks and Tasks matrix library [Parallel Comput. 57, 87 (2016)]. The expansion is based on second-order polynomials and accelerated by the scale-and-fold technique [J. Chem. Theory Comput. 7, 1233 (2011)]. We evaluate the performance of the implementation by computing the density matrix from the Fock matrix in the large-scale self-consistent field calculations. We demonstrate that the amount of communicated data per worker process tends to a constant with increasing system size and number of computer nodes such that the amount of work per worker process is fixed. (read more)