Multiplying two different density matrices, what physical situation will such need arises? This question arises from a reading in quantum chemistry:
In this link, the natural bonding orbitals (NBO) $\Theta_k$ (basically localised versions of molecular orbitals) are the eigenfunctions of a one electron reduced density matrix of an $n$ electron system $\Gamma (x,x') =n\int\cdots\int\psi^*(x,x_2,\cdots,x_n)\psi(x',x_2,\cdots,x_n)dx_2dx_3\cdots dx_n$, i.e.
\begin{equation}
\Gamma \Theta_k=p_k\Theta_k\tag{1}
\end{equation}
where $p_k$ are the occupation of each NBO thus $\sum_kp_k=n$
Now, imagine dividing $(1)$ both sides by $n$ we get $\frac{1}{n}\Gamma\Theta_k=P_k\Theta_k$, where $\sum_kP_k=1$
Translate all of this into bra-ket notation, it suggests something like this
$$\rho\lvert\psi\rangle=P_k\lvert\psi\rangle$$
where $\lvert\psi\rangle$ is a state vector and $\rho$ is a density matrix.
However, this is a bit strange as we never encountered any situation where both density matrices and state vector appearing at the same time, and that state vectors can always be written as a projector $\lvert\psi\rangle\langle\psi\rvert$ which is a density matrix.

1a. Is $\rho\lvert\psi\rangle$ a legal operation outside the context of eigenvalue equations?
1b. If it is legal, what physical context will it arises. how to interpret the resulting state vector?

Furthermore, using the above observations, we can also generalise it to something like this, written here in terms of matrix elements under some basis where $i,j,k$ is taken from said basis:
$\hat{\rho_1}\cdot\hat{\rho_2}=\sum_{j}\langle i\rvert \rho_1 \lvert j\rangle\langle j\rvert \rho_2 \lvert k\rangle$
where $\rho_1,\rho_2$ are two different density matrices
Now we knew that the purity of the mixed state is given by $\textrm{Tr}(\rho^2)$, thus we could ask whether there is something like $\textrm{Tr}(\rho_1\rho_2)$ or even $\textrm{Tr}(\rho_1 \ln \rho_2)$ and $\frac{\rho_1\rho_2\rho_1}{\textrm{Tr}(\rho_1\rho_2)}$
Since any density matrices can be written as a sum of projectors weighted by probability, the matrix multiplication $\rho_1\rho_2$ is effectively projecting one mixed state into another

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?
2b. Is $\rho_1\rho_2$, the multiplication of one density matrix onto a different one legal?
2c. What are examples of physical or experimental scenario will the need to project a mixed state to another mixed state arises?

 A: 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.
