$$Q=\int_{\text{all space}} \frac{\hbar \nu_g \mathbf{k}\mathbf{k}}{\exp[(\hbar \nu_g |\mathbf{k}|-\mathbf{k}\cdot\mathbf{u})/k_B T]-1}d\mathbf{k} $$
Please help me to integrate the above tensor expression in the infinite domain of the vector $k$.
I have tried to let $u$ in the direction of $k_z$ and then transform the current integral into a spherical coordinate with the following relation:
$$ k_x = k \sin \varphi \cos \theta $$ $$ k_y = k \sin \varphi \sin \theta $$ $$ k_z = k \cos \varphi $$
And I find $Q$ is a symmetric tensor, and $Q_{xy}=0,Q_{xz}=0,Q_{yz}=0$. However,$Q_{xx}=Q_{yy}\neq 0,Q_{zz}\neq 0$,and $Q_{xx}\neq Q_{zz}$. I guess $Q$ has relevance to the tensor $uu'$, the assumption that $u$ is in the direction of $k_z$ is just a simple case.
(This will not be a complete result so much as a strategy, since I don't recall all the ins-and-outs of this calculation.) While one can integrate directly, the more 'physical' strategy is to take advantage of the symmetry and tensorial form of the integral. First, we note that the integral is invariant under rotations about the origin; for convenience, then, we can orient our coordinate system so that $\mathbf{u}=u\hat{z}$ i.e. taking $\mathbf{u}$ to define the $z$-axis.
Second, we note that $Q$ is a rank-two spherically symmetric tensor. But the only such tensors available in the problem are the identity matrix $I$ and the outer product $\mathbf{u}\mathbf{u}=u^2\hat{z}\hat{z}$, so $Q$ must be a linear combination of these tensors. To this end it will be useful to introduce the longitudinal and transversal projection matrices $P_L:=\hat{u}\hat{u}=\hat{z}\hat{z}$ and $P_T:=I-P_L$ respectively. (Note that $P_L P_T=P_T P_L=0$). Then we assume that the $k$-space integral is of the form $Q(u) = Q_L(u)P_L+Q_T(u)P_T$ where $Q_L,Q_T$ are functions of $u$ which we need to find.
This structure can be exploited by taking traces. Since $\text{tr}(P_L)=1$ and $\text{tr} (P_T)=2$, we have
$$\text{tr}(P_L Q) = Q_L(u),\quad \text{tr}(P_T Q) = 2Q_T(u),$$
$$\implies Q(u) = \text{tr}(P_LQ)P_L+\frac{1}{2}\text{tr}(P_TQ)P_T$$
The problem is thereby reduced to evaluating the traces
$$\text{tr}(P_L Q)=\int_{\text{all space}} \frac{\hbar \nu_g k_z^2}{\exp[(\hbar \nu_g |\mathbf{k}|-k_z u)/k_B T]-1}d\mathbf{k},\\ \text{tr}(P_T Q)=\int_{\text{all space}} \frac{\hbar \nu_g (k_x^2+k_y^2)}{\exp[(\hbar \nu_g |\mathbf{k}|-k_z u)/k_B T]-1}d\mathbf{k}.$$ We have therefore traded one tensorial integral for two scalar integrals, a task best left for spherical coordinates. For now, I'll leave the remainder of the calculation to the reader.