  We calculate three-loop photon spectral density in QED with $N$ different
species of electrons. The obtained results were expressed in terms of iterated
integrals, which are either reduce to Goncharov's polylogarithms or can be
written in terms of one-fold integrals of harmonic polylogarithms and complete
elliptic integrals. In addition we provide threshold and high-energy
asymptotics of the calculated spectral density. It is shown, that the use of
the obtained spectral density correctly reproduces separately calculated
moments of corresponding photon polarization operator.
