Do the states in a decomposition $\rho=\sum_i p_i |\phi_i\rangle\!\langle \phi_i|$ need to be orthonormal? On Wikipedia it says:

Let $\mathcal H_S$ be a finite-dimensional Hilbert space, and consider a generic (possibly mixed) quantum state $\rho$ defined on $\mathcal H_S$, and admitting a decomposition of the form $\rho=\sum_i p_i|\phi_i\rangle \langle\phi_i|,$ for a collection of (not necessarily mutually orthogonal) states $|\phi_i\rangle\in\mathcal H_S$, and coefficients $p_i\ge 0 $ such that $\sum_i p_i=1$. Note that any quantum state can be written in such a way for some $\{|\phi_i\rangle\}$ and $\{p_i\}$.


*

*What is the propery of this set $\{|\phi_i\rangle\}$?


*

*Should it be normalized  $\langle\phi_i|\phi_i\rangle=1$?

*Should it be orthogonal  $\langle\phi_i|\phi_j\rangle=0$ for $i \neq j$?

*Should it ortho-normal $\langle\phi_i|\phi_j\rangle=\delta_{ij}$?



*Why does Wiki says that $\{|\phi_i\rangle\}$ not necessarily mutually orthogonal? If so, how can we assure that $\sum_i p_i=1$ in the non-orthogonal basis?

 A: You seem to be confusing the definition for a density matrix and the definition for an orthonormal basis. An orthonormal basis has to have the properties you described. A density matrix is something else. The terms $p_i$ can be interpreted as the probability to have the state $|\psi_i\rangle$, and this state can be arbitrary.
An example for a 2 level system. $p_1 = p_2 = 0.5$ and $|\psi_1\rangle = (1,0)$, and $|\psi_2\rangle = (1,1)/\sqrt 2$. This gives the perfectly valid density matrix $\rho = ((3/4, 1/4), (1/4, 1/4))$, even though the two states are not orthogonal
A: Remark: In order to avoid unnecessary mathematical complications, a finite dimensional Hilbert space (dimension $N$) will be assumed in the following discussion. The generalization to the infinite dimensional case is (in principle) straightforward.
Any convex linear combination $$\rho= \sum\limits_i p_i |\phi_i \rangle \langle \phi_i |, \quad p_i\ge 0, \, \sum\limits_i p_i=1 \qquad (1)$$
of an arbitrary set of pure states $| \phi_i\rangle$ (with $\langle \phi_i |\phi_i\rangle=1$) represented by the one-dimensional orthogonal projection operators $P_i= | \phi_i \rangle \langle \phi_i |$ describes the density operator $\rho$ of a possible mixed state. Note that the unit vectors vectors of the set $\{ |\phi_i \rangle \}_i$ do not have to be linearly independent nor do they have to form an orthonormal system. Even a convex linear combination in the form of an integral,
$$\rho =\int\limits \! d\alpha \,f(\alpha) |\phi_\alpha \rangle \langle \phi_\alpha|, \qquad f(\alpha) \ge 0, \ \int\! d\alpha \, f(\alpha) =1, \, \langle \phi_\alpha |\phi_\alpha \rangle =1 \qquad (2)$$ with a continuous probability distribution $f(\alpha)$ describes a possible density operator.
It is obvious that the representation $(1)$ of a given density operator $\rho$ is not unique. It is therefore convenient to find a certain "standard representation" for the density operator. As can be seen from $(1)$, a density operator is a positive selfadjoint operator ($\rho= \rho^\dagger \ge 0$) normalized as ${\rm Tr} \, \rho =1$. (Alternatively, these properties could be taken as the definition of a density operator.) As a consequence, the spectral theorem guarantees the existence of a complete orthonormal basis $\{ | \psi_n \rangle \}_{n=1}^N$ of eigenvectors, $\rho |\psi_n \rangle = \rho_n |\psi_n \rangle$ with nonnegative eigenvalues $\rho_n$ satisfying $\sum\limits_{n=1}^N \rho_n=1$. The spectral representation of the density operator is now given by $$\rho = \sum\limits_{n=1}^N \rho_n | \psi_n \rangle \langle \psi_n |, \quad \rho_n \ge 0, \, \langle \psi_n | \psi_m \rangle = \delta_{nm}, \, \sum\limits_{n=1}^N |\psi_n \rangle \langle \psi_n | = {\mathbf 1}_N. \qquad (3)   $$ Note that the representation chosen for the density operator $\rho$ of a given mixed state does, of course, not affect the physics. The expectation value of an oberservable $A$ in the state described by the density matrix $\rho$, given by $\langle A \rangle_\rho = {\rm Tr} (\rho A)$, is independent of the chosen decomposition of $\rho$.
P.S.: It is a nice homework exercise to find the spectral representation of the density matrix in the example given in the answer by peep!
A: It feels like point (2) in your question is the one that really confuses you. In a sense, you seem to think that $\operatorname{tr}\rho=1$ demands the states  $|\psi_i\rangle$ to form an orthonormal basis. Perhaps then the best way to approach this is to actually prove that $\operatorname{tr}\rho = \sum_i p_i$ regardless of what are the individual states in the decomposition.
Let therefore $\{|\psi_i\rangle\}$ be an arbitrary set of vectors, which we assume to be normalized (recall that in QM states are unit rays in a Hilbert space, so we should always normalize our kets). Define $$\rho=\sum_{i}p_i|\psi_i\rangle\langle \psi_i|$$
Let now $|e_n\rangle$ be an arbitrary orthonormal basis of the Hilbert space. The trace of $\rho$ is $$\operatorname{tr}\rho=\sum_n \langle e_n|\rho |e_n\rangle=\sum_{n,i}p_i\langle e_n|\psi_i\rangle \langle \psi_i|e_n\rangle$$
Now observe that rearranging and summing over $n$ first this is:
$$\operatorname{tr}\rho=\sum_i p_i \langle \psi_i|\left(\sum_n |e_n\rangle \langle e_n| \right)|\psi_i\rangle.$$
At this point we use the completeness of the  orthogonal basis and we observe that
$$\operatorname{tr}\rho=\sum_i p_i \langle \psi_i|\psi_i\rangle = \sum_i p_i,$$
where the last equality follows from the fact that the $|\psi_i\rangle$ are normalized. As such $\operatorname{tr}\rho=\sum_i p_i$ without the need for the $\{|\psi_i\rangle\}$ to form an orthonormal basis.
Finally, the states themselves cannot assure $\sum_i p_i=1$. You can pick any set of (normalized) kets you like and construct any density matrix you like by assembling $\rho = \sum_i p_i |\psi_i\rangle\langle \psi_i|$ where $\{p_i\}$ are any numbers restricted just to obey $\sum_i p_i=1$.
If the kets were initially not normalized, you just need to normalize them before constructing the mixed state.
