The main problem with the last statement,
$$
\begin{align}
C^{MP2} &= c\cdot c^T \\
\Rightarrow P^{MP2}_{ij} &= \sum_{a}C^{MP2}_{ia}(C^{MP2}_{ja}),
\end{align}
$$
is that one never has "MP2 MO coefficients" in this way; you would have MP2 natural orbital (NO) coefficients that come from diagonalizing this density, so it's the other way around:
$$
P^{\text{MP2}}C^{\text{MP2 NO}} = n^{\text{MP2 NO}}C^{\text{MP2 NO}},
$$
where $n^{\text{MP2 NO}}$ is the diagonal matrix of natural orbital occupation numbers. This still leaves the question of how to calculate $P^{\text{MP2}}$. Conceptually, the way I view it is that the set of $\{c_{ijab}^{(1)}\}$ are the amplitudes describing how to correct the uncorrelated ground-state wavefunction; they don't specify the wavefunction on their own, so $C^{\text{SCF}}$ will probably have to appear in the total density.
All of the following assumed a restricted wavefunction. $i,j,k,\dots$ are occupied MO indices, $a,b,c,\dots$ are virtual/unoccupied MO indices, $p,q,r,\dots$ are general (any) MO indices, and $\mu,\nu,\dots$ are general AO indices. The 2nd-order correction to the density, specifically the occupied-occupied and virtual-virtual blocks are given by Eqn. (11a) in Trucks [1]:
$$
\begin{align}
P_{ij}^{\text{MP2}} &= -\frac{1}{2} \sum_{kab} c_{jk}^{ab} c_{ik}^{ab} \\
P_{ab}^{\text{MP2}} &= +\frac{1}{2} \sum_{ijc} c_{ij}^{ac} c_{ij}^{bc},
\end{align}
$$
where I switched to using your notation. This is indeed similar to your idea from above; accumulate over both indices of the opposite space, then do a matrix multiplication.
The total density is then [Eqn. (7)]
$$
\mathscr{P} = P^{\text{SCF}} + P^{\text{MP2}},
$$
where $P_{pq}^{\text{SCF}} = \delta_{pq}$ because so far we're working in the MO basis. The back-transformation to the AO basis will be the same as for any other 2-index MO-basis quantity:
$$
\begin{align}
P_{\mu\nu} &= \sum_{pq}^{\text{all MOs}} C_{\mu p} P_{pq} (C_{\nu q})^{T} \\
\Rightarrow \mathbf{P}^{\text{AO}} &= \mathbf{C} \cdot \mathbf{P}^{\text{MO}} \cdot \mathbf{C}^{T}.
\end{align}
$$
Then, we can evaluate molecular properties that are represented as energy derivatives [Eqn. (6)]:
$$
\left.\frac{\partial \mathscr{E}}{\partial \lambda}\right|_{\lambda=0} = \sum_{\mu\nu} \mathscr{P}_{\mu\nu} Q_{\mu\nu}^{\lambda}
$$
where $\mathscr{E} = E_0^{(0)} + E_0^{(2)}$ and $Q_{\mu\nu}^{\lambda}$ is the matrix representation of the external perturbation $\lambda$ in the AO basis. This is analogous to the SCF result:
$$
\left.\frac{\partial E^{\text{SCF}}}{\partial \lambda}\right|_{\lambda=0} = \sum_{\mu\nu} P_{\mu\nu}^{\text{SCF}} Q_{\mu\nu}^{\lambda}
$$
Unrelaxed vs. relaxed density
The above definition of the density is incomplete, since it describes the unrelaxed correlated density. Our usual formulation of correlated methods is to describe correlation as some to-be-determined mixture of the virtual orbitals with the occupied orbitals. There is no such mixing in the above equations. In reality, $P_{ia}^{\text{MP2}} \overset{!}{\neq} 0$, so there is something that looks like a "transition density" with non-zero occupied-virtual blocks, in order to correct for the fact that our $C^{\text{SCF}}$ has no correlation. Allowing this mixing (which corresponds to orbital response) leads to the relaxed density. However, one has to solve the coupled-perturbed Hartree-Fock (CPHF or CPSCF) equations.
See between Eqns. (8) and (9) in [1] for the CPHF equations and how to solve them for the occupied-virtual mixing coefficients and finally the $P_{ia}^{\text{MP2}}$ block. Generally, note the similarities to the CIS density equations.
CPHF equations and first-order properties
As stated, when differentiating the energy with respect to some perturbation $\lambda$ to calculate first-order properties, one needs to contract the density with the operator property matrix. If you transform Eqn. (6) from the AO into the MO basis,
$$
\begin{align}
\left.\frac{\partial \mathscr{E}}{\partial \lambda}\right|_{\lambda=0} &= \sum_{pq} \mathscr{P}_{pq} Q_{pq}^{\lambda} \\
&= \sum_{pq} P_{pq}^{\text{HF}} Q_{pq}^{\lambda} + \sum_{pq} P_{pq}^{\text{MP2}} Q_{pq}^{\lambda},
\end{align}
$$
the SCF contribution simplifies due to the nature of the SCF density matrix in the MO basis,
$$
P_{pq}^{\text{HF}} = \left\{
\begin{array}{lr}
2 & \text{if }p = q\text{, p is an occupied MO} \\
0 & \text{otherwise}
\end{array}
\right.,
$$
leading to the expression
$$
\begin{align}
\left<\hat{\lambda}\right> &= \sum_{pq} P_{pq}^{\text{HF}} Q_{pq}^{\lambda} \\
&= 2 \sum_{i} Q_{ii}^{\lambda} \\
&= \mathrm{tr}\left(\mathbf{P}^{\text{AO}}\mathbf{Q}^{\text{AO}}\right) \\
&= \mathrm{tr}\left(\mathbf{P}^{\text{MO}}\mathbf{Q}^{\text{MO}}\right)
\end{align}
$$
for one-electron operators and wavefunctions that obey the Hellmann-Feynman theorem. I absorbed the 2 into $\mathbf{Q}^{\text{MO}}$, but this depends on the program implementations and their definition of the density. The point is that this contribution is relatively simple, and illustrates how for one-electron operators (dipole, quadrupole, angular momentum, ...), evaluating them for a converged Hellmann-Feynman wavefunction is very fast: one the density and property integrals are in the same basis, either do a matrix multiply and then trace, or do an element-wise (Schur) product and accumulate all the elements, so it scales as $O(N^3)$. No solution of CPHF equations is required. The same goes for the gradient with respect to nuclear coordinates $\frac{\partial E^{\text{SCF}}}{\partial x}$; it isn't a one-electron operator, and requires a mathematical trick to make $\frac{\partial P^{\text{SCF}}}{\partial x}$ disappear, but the SCF gradient doesn't require solving the CPHF equations either. This isn't true for $\frac{\partial E^{\text{MP2}}}{\partial x}$, which requires $\frac{\partial P^{\text{MP2}}}{\partial x}$, where you need to solve the CPHF equations to even form the full $P^{\text{MP2}}$.
Consider the remaining MP2 part of the property evaluation; the total density contains contributions from the occupied-occupied (oo) block, the virtual-virtual (vv) block, and the occupied-virtual (ov) block. I'm going to switch from $P^{\text{MP2}}$ and $c$ to $D$ and $t$ to make typing easier.
$$
\sum_{pq} D_{pq} Q_{pq}^{\lambda} = \sum_{ij} D_{ij} Q_{ij}^{\lambda} + \sum_{ab} D_{ab} Q_{ab}^{\lambda} + 2\sum_{ia} D_{ia} Q_{ia}^{\lambda}
$$
The first two terms can already be calculated from the unrelaxed density; $D_{ia}$ is the term that must come from solving the CPHF equations.
In general, the form of solving CPHF equations looks like $\mathbf{AU} = -\mathbf{V}$, where $\mathbf{A}$ is the orbital Hessian (2nd derivative of the energy with respect to orbital rotations), $\mathbf{V}$ is the gradient of some perturbation $\hat{V}$. In [1], it's $\lambda$, but Neese [3,4] and the Scandinavians usually use $V$ or $Q$ when the perturbation is associated with an observable. This perturbation gradient, termed the right-hand side (or RHS), describes how the perturbation might cause occupation of virtual orbitals, and the orbital Hessian describes how the response of one MO to some perturbation is coupled to the other MOs [5,6].
Start by rewriting $\sum_{ia} D_{ia} Q_{ia}^{\lambda}$ as the combination of a perturbation-dependent response $\mathbf{U}$ and the perturbation gradient $\mathbf{X}$:
$$
2\sum_{ia} D_{ia} Q_{ia}^{\lambda} = 2\sum_{ia} X_{ai} U_{ai}^{\lambda}
$$
$\mathbf{X}$ serves the same purpose as the $\mathbf{V}$ mentioned above, but rather than coming from an external perturbation, it is an intermediate quantity that describes the MP2 wavefunction.
First, the perturbation-dependent response. This part is not specific to MP2, and is also important for higher-order SCF-based molecular response [3].
$$
\begin{align}
U_{ai}^{\lambda} &= \sum_{bj} (\mathbf{A}^{-1})_{ai,bj} Q_{bj}^{\lambda} (\varepsilon_{j} - \varepsilon_{b})^{-1} \\
A_{ai,bj} &= 1 + \left( \left<ab||ij\right> - \left<aj||ib\right> \right) (\varepsilon_{a} - \varepsilon_{b})^{-1}
\end{align}
$$
The $\varepsilon$ are MO energies. If we don't assume the use of canonical MOs, this energy denominator is actually $F_{ab}S_{ij} - F_{ij}S_{ab}$.
Jumping ahead, the relationship between the ov density, the orbital Hessian, and the MP2 RHS is [Eqn. (9)]
$$
\sum_{bj} D_{bj} (\varepsilon_{j} - \varepsilon_{b}) A_{bj,ai} = X_{ai}
$$
See how this is structurally similar to solving for $U$, but the RHS is different. To show the equivalence between all quantities,
$$
\begin{align}
2\sum_{ai} X_{ai} U_{ai}^{\lambda} &= 2\sum_{ai} X_{ai} \sum_{bj} (\mathbf{A}^{-1})_{ai,bj} Q_{bj}^{\lambda} (\varepsilon_{j} - \varepsilon_{b})^{-1} \\
&= 2\sum_{bj} D_{bj} Q_{bj}^{\lambda}.
\end{align}
$$
The explicit equation for $X_{ai}$ is [Eqn. (11b)]
$$
X_{ai} = \sum_{jk} D_{kj} \left<ki||ja\right> + \sum_{bc} D_{bc} \left<bi||ca\right> - \frac{1}{2} \left[ t_{kj}^{ab} \left<ib||kj\right> + t_{ij}^{cb} \left<cb||aj\right> \right],
$$
which is in fact the MP2 Lagrangian, $L_{ai}$ [2, Eqn. (11)]. Again, notice the similarity to the CIS Lagrangian. The Lagrangian is required because the wavefunction doesn't need to be relaxed in order to calculate the correlation energy, so you pay the price here. These equations must be evaluated for MP2 geometry optimizations (during the gradient part) or for fully-relaxed MP2 molecular properties. To be clear, if one is only interested in the density and not molecular properties, only $\mathbf{DA} = \mathbf{X}$ needs to be solved for, and not the equations involving $\mathbf{U}$. In terms of scaling, SCF-based CPHF is about $O(N^4)$ (after back-transforming the $\mathbf{A}^{-1}\mathbf{Q}$ terms to be AO-direct) but might require fewer iterations than the SCF itself. According to [2], the MP2 gradient is something like $O(ON^4) + O(O^{2}V^{3})$, but they can do similar tricks during iterations to cut down on the number of fifth-order steps. In my experience, MP2 gradient calculations are about as expensive as an SCF polarizability calculation.
References
Trucks, G. W.; Salter, E. A.; Sosa, C.; Bartlett, R. J. Theory and implementation of the MBPT density matrix. An application to one-electron properties. Chem. Phys. Lett. 1988, 147 (4), 359-366. DOI: 10.1016/0009-2614(88)80249-5
Frisch, M. J.; Head-Gordon, M.; Pople, J. A. A Direct MP2 gradient method. Chem. Phys. Lett. 1990, 166 (3), 275-280. DOI: 10.1016/0009-2614(90)80029-D
Neese, F. Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling. Coordination Chemistry Reviews 2009, 253 (5-6), 526-563. DOI: 10.1016/j.ccr.2008.05.014
I cannot recommend the above paper enough; it contains excellent descriptions of how electric and magnetic response are calculated with explicit equations that are actually usable. There, you can see how CPHF is used for molecular response properties and the connection with excitation energies, found by solving $\mathbf{AU} = \mathbf{0}$ (diagonalization rather than solving a set of linear equations).
Toulouse, J. Introduction to the calculation of molecular properties by response theory.
Parrish, R. M. Psi4 Development: Theory Supplement, CPHF for Dipole Polarizabilities.
