If $|\psi\rangle$ is a n-electron state and
$$
\psi(x_1, x_2, x_3, \dots, x_n) = \langle x_1, x_2, x_3, \dots, x_n | \psi\rangle
$$
is its (anti-symmetric) n-electron wave function in position representation, then
$$
\frac{1}{n} {\hat\Gamma} = Tr_{2,3,\dots, n} |\psi\rangle \langle \psi| = \int{dx_1 dx'_1\,|x_1\rangle \,\Big( \int{dx_2 dx_3 \dots dx_n \psi^*(x_1, x_2, x_3, \dots, x_n) \psi(x'_1, x_2, x_3, \dots, x_n)\; }\Big) \langle x'_1| }
$$
is the single-electron reduced density matrix and
$$
\frac{1}{n}\Gamma(x, x') = \langle x | {\hat \Gamma} | x' \rangle
$$
gives its matrix elements in position representation. Notice that since $\psi(x_1, x_2, x_3, \dots, x_n)$ is anti-symmetric in the electron coordinates, it doesn't matter which electrons are traced out.
The equation
$$
\frac{1}{n} {\hat\Gamma} |\Phi_k\rangle = p_k |\Phi_k\rangle
$$
is indeed the eigenvalue equation for the reduced single-electron density matrix ${\hat \Gamma}$ and the NBOs $|\Phi_k\rangle$ are its eigenstates. This means that an electron can be found with probability $p_k$ in $|\Phi_k\rangle$. Or in density functional theory (DFT) lingo, the (possibly fractional) occupation number of orbital $\Phi_k$ is $n_k = np_k$. Otherwise, everything that pertains to (reduced) density matrices applies to ${\hat \Gamma}$ as well. So,
1a. Is $\rho|\psi\rangle$ a legal operation outside the context of eigenvalue equations?
Yes, absolutely, $\rho$ is just another operator.
1b. If it is legal, what physical context will it arises. how to interpret the resulting state vector?
It doesn't necessarily mean that $\rho|\psi\rangle$ has some special physical significance. On the other hand, for any state $|\psi\rangle$ the matrix element $\langle \psi |\rho |\psi\rangle$ has the meaning of "probability to measure the system (electron) in the pure state $|\psi\rangle$ when it is in the mixed state $\rho$". For ${\hat \Gamma}$ this means $\langle \Phi |\frac{1}{n}{\hat \Gamma} |\Phi\rangle$ gives the probability to find an electron in orbital $|\Phi\rangle$ if the reduced single-electron state is $\frac{1}{n}{\hat \Gamma}$.
2a. Since a density matrix is a matrix representation of the density operator, given any linear operator ${\hat O}$ that is used in quantum mechanics and is not necessary an observable nor the unitary time evolution operator $U(t,t′)$, is the product ${\hat O}\rho$ legal? What is its meaning?
Yes, the product ${\hat O}\rho$ is well-defined, but again, it doesn't have special physical significance by itself.
2b. Is $\rho_1\rho_2$, the multiplication of one density matrix onto a different one legal?
Again, yes. But it usually doesn't arise when $\rho_1$ and $\rho_2$ both refer to the same single electron coordinates (it is tempting to say "the same electron", but then electrons are indistinguishable, so we can't possibly know that it is "the same"). It usually means a direct product $\rho_1 \otimes \rho_2$ for two electrons (two sets of electron coordinates).
2c. What are examples of physical or experimental scenario where the need to project a mixed state to another mixed state arises?
It's not a projection!
Technically all quantum mechanical equations are linear in the density matrix of the system, so whenever you see this sort of product you know two things:
there are (at least) two electrons involved and the (direct) product $\rho_1 \otimes \rho_2 $ refers to (part of) a two-electron state ;
very likely there is some self-consistent approximation about the nature of the two-electron reduced density matrix (obtained by tracing the n-electron density matrix $|\psi\rangle\langle \psi|$ over all but two electrons). This is because the two-electron density matrix is a product of single electron density matrices only when the electrons are uncorrelated.
In DFT matrix elements of products for $\rho_1 = \rho_2 = \rho$ are ubiquitous and you can easily find them in a variety of quantities: the Hartree energy, the exchange energy, the two-electron correlation function, the exchange correlation function, etc.
Best Answer
Proposition: Let $\{ \left| \psi ^{(k)}\right> :k\}$ be any collection of normalized vectors in a hilbert space, let $0\leq p_k\leq 1$ be such that $\sum _kp_k=1$, and define $\rho :=\sum _kp_k\left| \psi ^{(k)}\right> \left< \psi ^{(k)}\right|$. Then, (i) $\operatorname{tr}(\rho )=1$ and (ii) $\operatorname{tr}(\rho ^2)\leq 1$.
Remark: "Any" really means any: this collection needs to be neither orthonormal nor a basis (though each $\left| \psi ^{(k)}\right>$ definitely needs to be normalized!).
Proof: Let $\{ \left| e_i\right> :i\}$ be an orthonormal basis of the hilbert space (not necessarily eigenvectors for $\rho$), and write $\left| \psi ^{(k)}\right> =\sum _ic^{(k)}_i\left| e_i\right>$. (Note that the index $i$ runs over (in general) very different index set that $k$ does.) Then, $$ \operatorname{tr}(\rho ):=\operatorname{tr}\left( \sum _kp_k\left| \psi ^{(k)}\right> \left< \psi ^{(k)}\right|\right) =\sum _kp_k\operatorname{tr}\left( \left| \psi ^{(k)}\right> \left< \psi ^{(k)}\right| \right)=\sum _kp_k\cdot 1=1, $$ where one can see that $\operatorname{tr}\left( \left| \psi ^{(k)}\right> \left< \psi ^{(k)}\right| \right) =1$ by picking any orthonormal basis of the hilbert space which contains $\left| \psi ^{(k)}\right>$ and using the definition of the trace.
Hence, $\rho$ is trace-class (Disclaimer: You actually need to check that $\operatorname{tr}(|\rho |)<\infty$, but it shouldn't be too hard to show that $|\rho |=\rho$.), hence compact, and so there is an orthonormal basis $\{ \left| f_i\right> :i\}$ consisting of eigenvectors of $\rho$: $\rho \left| f_i\right> =\lambda _i$. Expanding the equation $\lambda _i=\left< f_i\right| \rho \left| f_i\right>$ shows that $0\leq \lambda _i\leq 1$, and so $$ \operatorname{tr}(\rho ^2)=\sum _i\left< f_i\right| \rho ^2\left| f_i\right> =\sum _i\lambda _i^2\leq \sum _i\lambda _i=1. $$ $\square$