Why is density fitting not used to reduce four-center integrals further down to two-center integrals?
I try to answer your question, although I'm not completely sure if I understand what you are asking for.

You want a decomposition like for example:
$$
 (pq|rs) = \sum_P A^p_P A^q_P A^r_P A^s_P
$$
am I correct? 

In the framework of density fitting this type of decomposition is not possible without further steps/tricks.
If we look at the two electron integrals
\begin{equation}
(pq|rs)=\int d\pmb{r}_1 \int d\pmb{r}_2 \phi_p(\pmb{r}_1) \phi_q(\pmb{r}_1) \frac{1}{\pmb{r}_{12}} \phi_r(\pmb{r}_2) \phi_s(\pmb{r}_2) \quad,
\end{equation}
we can consider this electron repulsion integral as repulsion between two generalized electron densities
\begin{equation}
(pq|rs)=\int d\pmb{r}_1 \int d\pmb{r}_2 \rho_{pq}(\pmb{r}_1)  \frac{1}{\pmb{r}_{12}} \rho_{rs}(\pmb{r}_2) \quad,
\end{equation}
with $\rho_{pq}(\pmb{r}) = \phi_p(\pmb{r}) \phi_q(\pmb{r})$ and $\rho_{rs}(\pmb{r}) = \phi_r(\pmb{r}) \phi_s(\pmb{r})$.
This density maybe approximated with the aid of an auxiliary basis set:
\begin{equation}
 \tilde{\rho}_{pq}(\pmb{r}) = \sum^{N_{aux}}_P d_P^{pq} \chi_P(\pmb{r})
\end{equation}
Hence the name "density fitting". The coefficients $d_P^{pq}$ are determined by minimizing the functional
\begin{equation}
  \Delta_{pq} = \int d\pmb{r}_1 \int d\pmb{r}_2 \frac{\left[ {\rho}_{pq}(\pmb{r_1}) - \tilde{\rho}_{pq}(\pmb{r_1}) \right] \left[ {\rho}_{rs}(\pmb{r_2}) - \tilde{\rho}_{rs}(\pmb{r_2}) \right]}{r_{12}} 
\end{equation}
and leads according to Dunlap [1] to:
\begin{equation}
 d_Q^{pq} = \sum_{P} \left(pq|P\right) \left[\pmb{J}^{-1}\right]_{PQ} \quad,
\end{equation}
with
\begin{equation}
 \left(pq|P\right) = \int d\pmb{r}_1 \int d\pmb{r}_2 \phi_p(\pmb{r}_1) \phi_q(\pmb{r}_1) \frac{1}{\pmb{r}_{12}} \chi_P(\pmb{r}_2)
\end{equation}
\begin{equation}
 J_{PQ} = \int d\pmb{r}_1 \int d\pmb{r}_2 \chi_P(\pmb{r}_1) \frac{1}{\pmb{r}_{12}} \chi_Q(\pmb{r}_2)
\end{equation}
Then we can re-write the electron repulsion integral as
\begin{equation}
 (pq|rs)=\sum_{PQ} \left(pq|P\right) \left[\pmb{J}^{-1}\right]_{PQ} \left(Q|rs\right)
\end{equation}
The framework of density fitting does not give us the possibility to decompose $\left(pq|P\right)$ further,
but it is possible in combination with more general tensor decomposition techniques like the canonical decomposition (CP)[2],
in which we approximate a tensor in the format you were fishing for (as long as I understood you correctly)
\begin{equation}
 X_{i,j,...d} = \sum_{r}^{R} \omega_r  a_{ir}^{(1)} a_{jr}^{(2)} ... a_{dr}^{(d)} \quad,
\end{equation}
where $R$ is the rank of the tensor. Without losing generality, we will in the following absorb the weighting factor $\omega_r$
in one of the factors.

We can apply this general tensor decomposition to the density fitting representation of the two electron integral.
For the $[J^{-1}]_{PQ}$ matrix we can simply us a SVD
\begin{equation}
 [J^{-1}]_{PQ}  \approx \sum_{r_x}^{R_x} \bar{U}^{P}_{r_x} \bar{V}^{Q}_{r_x}  \quad,
\end{equation}
where the singular values $\sigma^{1/2}$ are absorbed into the vectors
For the 3-index integrals we do a CP decomposition
\begin{equation}
  \left(pq|P\right) \approx \sum_{r}^{R} A_{r}^{p} A_{r}^{q} A_{r}^{P}
\end{equation}
Thus the 4-index integrals can be written as
\begin{equation}
 \left(pq|rs\right)  = \sum_{PQ} \left( \sum_{r_1}^{R_1} A_{r_1}^{p} A_{r_1}^{q} A_{r_1}^{P} \right)
                                 \left( \sum_{r_x}^{R_x} \bar{U}^{P}_{r_x} \bar{V}^{Q}_{r_x} \right)
                                 \left( \sum_{r_2}^{R_2} A_{r_2}^{r} A_{r_2}^{s} A_{r_2}^{Q} \right) \quad.
\end{equation}
In some sense we are done since only two index quantities occur, but we can simplify the expressions a bit.
The sums can be rearranged to
\begin{equation}
  \left(pq|rs\right) = \sum_{r_1}^{R_1}\sum_{r_2}^{R_2}\sum_{r_x}^{R_x} \left( \sum_{P} A_{r_1}^{P}  \bar{U}_{r_x}^P  \right) \left( \sum_{Q} \bar{V}_{r_x}^Q   A_{r_2}^{Q} \right)
                         A_{r_1}^{p} A_{r_1}^{q} A_{r_2}^{r} A_{r_2}^{s} \quad,
\end{equation}
and a summation over $P$ and $Q$ yields
\begin{equation}
  \label{eqn:thc1}
  \left(pq|rs\right) = \sum_{r_1}^{R_1}\sum_{r_2}^{R_2}\sum_{r_x}^{R_x} P_{r_1 r_x} Q_{r_x r_2}
                         A_{r_1}^{p} A_{r_1}^{q} A_{r_2}^{r}  A_{r_2}^{s} \quad.
\end{equation}
Only $P_{r_1 r_x}$ and $Q_{r_x r_2}$ depend on $r_x$ and via summation we introduce a new coefficient matrix $B_{r_1 r_2}$, which
combines the expansions length $R_1$ and $R_2$.
\begin{equation}
  \left(pq|rs\right) = \sum_{r_1}^{R_1}\sum_{r_2}^{R_2} B_{r_1 r_2} A_{r_1}^{p} A_{r_1}^{q} A_{r_2}^{r} A_{r_2}^{s} \quad,
\end{equation}
which is quite similar to the THC format introduced by Martinez and co-workers.[3]
In one additional step we obtain a CP factorization for the full electron repulsion integral
\begin{equation}
  \left(pq|rs\right) = \sum_{P}^{R_1\cdot R_2} \bar{A}_{P}^{p}  \bar{A}_{P}^{q} \bar{A}_{P}^{r} \bar{A}_{P}^{s} 
\end{equation}
So finally we arrived at a decomposition in only two index quantities and there are several groups[4,3,5], who looked
at this integral format(s). It is not of widespread use since so far the available algorithms for the required CP decomposition
are quite costly. In combination with density fitting the costs will scale with $\mathcal O(N^4)$, which is not that bad,
but the prefactor is so large that for a long range of molecular sizes a CCSD(T) calculation will finish quicker than the decomposition
of the two electron integral tensor.

The story is a bit different for the intermediate step, which can be identified as the THC format. Martinez and co-workers
are heavily working on using this format and enabled $\mathcal O(N^4)$ scaling implementations for MP2 and CCSD. However,
for MP2 it is to the best of my knowledge not really compatible with DF-MP2 implementations.

So after this long answer the short conclusion: Such integral formats exist in form of the THC format or CP decomposed integrals,
but further work is necessary to increase the efficiency and resolve some other problems.



References:


Dunlap, B. I., Connolly, J. W. D. &amp; Sabin, J. R. On Some Approximations in Applications of X Alpha Theory. The Journal of Chemical Physics 71, 3396–3402 (1979). 
Carroll, J. D. &amp; Chang, J.-J.   Analysis of individual differences in multi-dimensional scaling via an n-way generalization of “eckartyoung” decomposition. Psychometrika 35, 283–319 (1970). URL http://dx.doi.org/10.1007/BF02310791
Hohenstein, E. G., Parrish, R. M. &amp; Martínez, T. J. Tensorhypercontraction density fitting. I. quartic scaling second- and third-order Møller-plesset perturbation theory. The Journal of Chemical Physics 137 (2012). URL http://scitation.aip.org/content/aip/journal/jcp/137/4/10.1063/1.4732310.
Benedikt, U., Auer, A. A., Espig, M. &amp; Hackbusch, W. Tensor decomposition in post-Hartree Fock methods. I. Two-electron integrals and MP2. The Journal of Chemical Physics 134 (2011). URL http://scitation.aip.org/content/aip/journal/jcp/134/5/10.1063/1.3514201.
Schmitz, G., Madsen, N. K. &amp; Christiansen, O. Atomic-batched tensor decomposed two-electron repulsion integrals. The Journal of Chemical Physics 146, 134112 (2017). URL https://doi.org/10.1063/1.4979571. https://doi.org/10.1063/1.4979571