How do we arrive that the form of the thermal density matrix? I'm following the notes by Pieter Kok, the relevant material for this question is around page 29. I learned about the density matrix for just a single spin-1/2 particle with unknown initial conditions, knowing only the  probabilities to be in the state in $\left|0\right>$ with classical probability 1/2 and in $\left|+\right>$ with classical probability 1/2, the density matrix is
$$\rho=\dfrac{1}{2}\left|0\right> \left<0\right|+\dfrac{1}{2}\left|+\right> \left<+\right|\tag{1}.$$
For more general cases one can define the density matrix as
$$\rho=\sum_kP_k\left|k\right>\left<k\right| \tag{2},$$
where $P_k$ is the classical probability to be in the state $\lvert k\rangle$.
But in quantum statistical mechanics we define the density matrix for the grand canonical ensemble as
$$\rho=\frac{1}{\mathrm{Tr}[\mathrm{e}^{-\beta(H-\mu N)}]}\mathrm{e}^{-\beta(H-\mu N)} \tag{3}.$$
How can I understand the relation between these two definitions? I am guessing this will be related to how we deal with the initial conditions for few-body systems and many-body systems.
Another part of my question comes from the derivation of linear response theory in the notes written by Sunkai, particular in section 3.2. For the thermal average of some observable $B$ one can use
$$\left<B\right>=Tr[\rho B] \tag{4},$$
where $\rho$ is the density matrix defined in eq. (3). Note that eq. (4) is expressed in the Schrödinger picture. One can transform this definition into the Heisenberg picture for the purpose of calculation: 
$$\left<B(t)\right>=Tr[\color{red}{\rho(t)} B(t)] \tag{5}.$$
For the red part I am very confused. Why one can still talk about the evolution of density matrix in Heisenberg picture, as done in the notes? Ìt should be constant because for the thermal mixed state the density matrix plays the same role just as the state vector. So am I misunderstanding something?
 A: The relation between the general definition of 
$$ \rho = \sum_k P_k \lvert k\rangle\langle k\rvert$$
as the density matrix for a mixed state with probabilities $P_k$ for the states $\lvert k\rangle$ and the thermal density matrix is straightforward:
Rather by definition, the grand canonical ensemble assigns classically to each microstate the probability
$$ P(E,N) = \frac{1}{Z}\mathrm{e}^{\beta(\mu N - E)}$$
and if we now take an eigenbasis $\lvert x_i\rangle$ of the self-adjoint operator $\mu N-H$, $P$ is simply a function of these $x_i$ and we have
$$ \rho = \sum_i P(x_i)\lvert x_i\rangle\langle x_i\rvert = \frac{1}{Z}\sum_i \mathrm{e}^{\beta x_i}\lvert x_i\rangle\langle x_i\rvert = \frac{1}{Z}\mathrm{e}^{\beta (\mu N - H)},$$
since every self-adjoint operator may be written as $O = \sum_k O_k \lvert o_k\rangle\langle o_k\rvert$ for $\lvert o_k\rangle$ an eigenbasis with eigenvalues $O_k$ - this is just the statement that these operators are diagonal in their eigenbasis.
As for your second question, the density matrix is not time-dependent in the Heisenberg picture. It is time-dependent in the Schrödinger picture and obeys a "wrong-signed" Heisenberg equation of motion, the von Neumann equation, see also this question.
