This procedure is analogously done for the nuclear two-particle 4-index K i 0 j 0 a 0 b 0 integrals. However, for the two-particle 4-index integrals built from electrons and nuclei K ii 0 aa 0 the density fitting can be employed within the electronic or nuclear auxiliary basis set. This was already further explored by the group of Hammes-Schi↵er. It was found that by employing the density fitting approximation with the nuclear density fitting basis the overall error in forming those integrals can be significantly reduced. Therefore, we also employ the nuclear density fitting within the nuclear part of the mixed interactions between electrons and nuclei. Moreover, we unify both the electronic and nuclear auxiliary fitting basis sets since this avoids the recomputation and repeated transformation of the 3-index integrals within the electronic part. This approach reduces the overall computational costs, especially if the electronic subsystem exceeds the nuclear subsystem, which we have discussed previously in the construction of our LDF-NEO-RHF program. Overall, the asymptotic scaling with the fifth order is not reduced by the density fitting approximations within the DF-NEO-MP2 program. However, the density fitting leads to a low prefactor of the nominally quintic scaling which leads to a significant computational speed-up for small to mid-sized molecules as we will demonstrate later.