1RDMs are just simple matrices, which can be stored in either dense or sparse form, possibly combined with triangular storage (the matrix is often symmetric). 1RDMs are passed around in a number of formats, like Gaussian formatted checkpoint and Molden, and can be visualized for density isosurfaces etc.
2RDMs are a bit more problematic, since they can become very large; because of this it might be that not all programs even store them in full but rather form and process only selected blocks at a time. The 2RDMs also contain more symmetries due to the particle interchanges, e.g. $\Gamma(pq|rs)=-\Gamma(qp|rs)=-\Gamma(pq|sr)=\Gamma(qp|sr)$ or $\Gamma(pq|rs)=\Gamma(rs|pq)$ which means that there may be several choices on how to do this.
However, the assumption here appears to be that the system is so small that all the $O(N^4)$ integrals can be stored. In this case, one might think about using the FCIDUMP format to store the density matrices, too. The original FCIDUMP format doesn't store the Fock matrix, but extended versions do; this means that you should be able to save the 1RDM using a Fock matrix convention and the 2RDM using the two-electron integral convention.