Is there a relation between the density matrix and the density of the position probability? Are they the same concept? We have a system of two bosons particles and we are interested in calculating the one-particle density and two-particle-density when both are in different states.
So, to do that, I consider the following:
First, we know that the exchange of any two identical bosons must be symmetric, therefore the that  wave function for a two-particle bosons is given by:
$$ \psi^s\left(\vec{r_1},\vec{r_2} \right) = \frac{1}{2}\left[ \psi_a\left(\vec{r_1}\right) \psi_b\left(\vec{r_2}\right) +\psi_a\left(\vec{r_2}\right) \psi_b\left(\vec{r_1}\right)   \right] $$
Here is my doubt:
If the system is composed of two bosons how can I calculate the one-particle density, Which density does it refer to?
... $\\$
For two- particles:
$$ \left| \psi^s\left(\vec{r_1},\vec{r_2} \right) \right|^2 $$
In many books, I found the concepts of density operator and density probability? Are they the same?
 A: Assuming your single particle states are orthonormal and $a$ and $b$
are different, the prefactor for $\psi^s$ should be $\frac{1}{\sqrt{2}}$
for a normalized state.
The density matrix is the general quantum form for the case where your
system can be in a mixed state. That is, it can be thought of
as describing an ensemble of systems, each in a pure quantum state with
a given probability. For a pure system, like you describe, you can
use either normal expectation values or traces over the density matrix
and the operator. The diagonal elements of the one-body reduced density
matrix in position space will be the one-body density.
Just using expectation values is the simplest way to get the
density for your problem.
The position probability density is the expectation value
of the one-particle density operator. For a system with $N$ particles
that would be
\begin{equation}
\hat \rho (\vec r) = \sum_{i=1}^N \delta^3(\vec r- \hat{\vec r_i})
\end{equation}
where $\vec r$ is the position you measure the density, and
$\hat{\vec r_i}$ are the position operators. That is, the density operator
when integrated over a volume, should count the number of particles in
the volume.
The one-particle density is
for your case the expectation value
\begin{equation}
\langle \hat \rho(\vec r)\rangle = \int d^3r_1 d^3r_2
\left [ \delta^3(\vec r-\vec r_1)+\delta^3(\vec r-\vec r_2)\right ]
|\psi^s(\vec r_1,\vec r_2)|^2
= |\psi_a(\vec r)|^2+|\psi_b(\vec r)|^2
\end{equation}
