By averaging over positions of all but one of the electrons, the one-electron density matrix for electronic eigenfunctions is a convenient object in understanding how electrons are distributed around an atom or molecule. We rigorously study the regularity of this density matrix. Firstly, we show that the density matrix is analytic away from the nuclear positions and the diagonal. This analyticity is not immediate since the underlying wavefunction has points of non-smoothness which are integrated over. However, using certain directional derivatives known as cluster derivatives we can differentiate along the singularities to avoid them contributing adversely to the integral. The Coulomb potential, which is analytic away from a set of singularities, has certain analyticity bounds. These can be inherited by the wavefunction solutions to the time-independent Schrödinger equation using elliptic regularity. These bounds are then conferred onto derivatives of the density matrix in order to prove its analyticity. The density matrix is then studied at the diagonal, where pointwise bounds are obtained to derivatives both in the direction along the diagonal (u-derivatives) and perpendicular to it (v-derivatives). We find that up to four v-derivatives of the density matrix may be taken for the function to remain bounded in the vicinity of the diagonal. Whereas arbitrarily many u-derivatives may be taken without contributing to a worsened singularity at the diagonal. To prove this result, we state and prove a new pointwise bound for cluster derivatives of wavefunctions where multiple clusters are involved. This bound separates the contributions from each cluster and is likely to be of independent interest.