We evaluate the density matrix of an arbitrary quantum mechanical system in terms of the quantities pertinent to the solution of the time-dependent density functional theory (TDDFT) problem. Our theory utilizes the adiabatic connection perturbation method of G\"orling and Levy, from which the expansion of the many-body density matrix in powers of the coupling constant $\ensuremath{\lambda}$ naturally arises. We then find the reduced density matrix ${\ensuremath{\rho}}_{\ensuremath{\lambda}}(\mathbf{r},{\mathbf{r}}^{\ensuremath{'}},t)$, which, by construction, has the $\ensuremath{\lambda}$ independent diagonal elements ${\ensuremath{\rho}}_{\ensuremath{\lambda}}(\mathbf{r},\mathbf{r},t)=n(\mathbf{r},t)$,
n(r,t) being the particle density. The off-diagonal elements of ${\ensuremath{\rho}}_{\ensuremath{\lambda}}(\mathbf{r},{\mathbf{r}}^{\ensuremath{'}},t)$ contribute importantly to the processes unaccessible via the density, directly or by the use of the known TDDFT functionals. Of those, we consider the momentum-resolved photoemission, doing this to the first order in $\ensuremath{\lambda}$, i.e., on the level of the exact exchange theory. In illustrative calculations of photoemission from the quasi-2D electron gas and isolated atoms, we find quantitatively strong and conceptually far-reaching differences with the independent-particle Fermi's golden rule formula.