Let us assume a linear mixed model of the form
$$ \boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{Z}\boldsymbol{b} + \boldsymbol{\epsilon} $$
where $\boldsymbol{X}$ is the fixed-effect design matrix, $\boldsymbol{\beta}$ the fixed-effect coefficient vector, $\boldsymbol{Z}$ a block-diagonal random-effect design matrix, $\boldsymbol{b}$ the corresponding random effects and $\boldsymbol{\epsilon}$ the idiosyncratic error vector. As usual, the idiosyncratic errors are assumed i.i.d. and the random effects are distributed according to covariance matrix $\boldsymbol{Q}$.
$$\boldsymbol{\varepsilon}_i\left|\boldsymbol{X}_i \sim \mathrm{N}\left(\mathbf{0}, \sigma^2 \boldsymbol{I}_{n_i}\right), \quad \boldsymbol{b}_i\right| \boldsymbol{X_i} \sim \mathrm{N}(\mathbf{0}, \boldsymbol{Q})$$
Based on these assumptions, we obtain the full idiosyncratic error matrix $\boldsymbol{Q}$ and the full random effect covariance matrix $\boldsymbol{D}$ and $\boldsymbol{V}$, which is the covariance matrix of $\boldsymbol{y}$.
$$\boldsymbol{R}=\operatorname{blockdiag}\left(\sigma^2 \boldsymbol{\Sigma}_{n_1}, \ldots, \sigma^2 \boldsymbol{\Sigma}_{n_i}, \ldots, \sigma^2 \boldsymbol{\Sigma}_{n_m}\right)$$
$$\boldsymbol{D}=\operatorname{blockdiag}(\boldsymbol{Q}, \ldots, \boldsymbol{Q}, \ldots, \boldsymbol{Q})$$
$$ \left(\begin{array}{l} \boldsymbol{b} \\ \boldsymbol{\varepsilon} \end{array}\right) \sim \mathrm{N}\left(\left(\begin{array}{l} \boldsymbol{0} \\ \boldsymbol{0} \end{array}\right),\left(\begin{array}{ll} \boldsymbol{D} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{R} \end{array}\right)\right) $$
$$\quad$$
$$\boldsymbol{V}=\boldsymbol{R}+\boldsymbol{Z} \boldsymbol{D} \boldsymbol{Z}^{\prime}$$
Unbiased estimates of the random effects can then be obtained according to
$$ \hat{\boldsymbol{b}} = \boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\hat{\beta}}) $$
$$\quad$$
I am interested in the covariance matrix of the estimated random-effect vector $\boldsymbol{b}$ and the unobservable idiosyncratic error vector $\boldsymbol{\epsilon}$, i.e. I want to obtain
$$ Cov[\hat{\boldsymbol{b}}, \boldsymbol{\epsilon}] $$
So far I tried the following
$$ Cov[\hat{\boldsymbol{b}}, \boldsymbol{\epsilon}] = E[(\boldsymbol{\hat{b}} - E[\boldsymbol{\hat{b}}])(\boldsymbol{\epsilon} - \boldsymbol{0})^{\prime}] = E\left[\left(\boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\hat{\beta}}) - \boldsymbol{b}\right)\boldsymbol{\epsilon}^{\prime}\right] = $$ $$ E\left[\left(\boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}\boldsymbol{y}- \boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}\boldsymbol{X}\boldsymbol{\hat{\beta}}\right)\boldsymbol{\epsilon}^{\prime} - \boldsymbol{b}\boldsymbol{\epsilon}^{\prime} \right] = E\left[\left(\boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}\boldsymbol{y}- \boldsymbol{\hat{D}}\boldsymbol{Z}^{\prime}\boldsymbol{V^{-1}}\boldsymbol{X}\boldsymbol{\hat{\beta}}\right)\boldsymbol{\epsilon}^{\prime}\right] $$
I do not know how to continue from here and I am even not quite sure if I got everything right until here.
Could anybody help me out? I would greatly appreciate your help!