Since the multimode Fock states $|\mathbf{n}\rangle=|n_1\rangle\otimes |n_2\rangle\otimes\cdots\otimes |n_M\rangle$ form an orthonormal basis, one can treat each state by their probability distribution alone, with elements
$p_{\mathbf{n}}=\prod_{j=1}^M P(n_j|\mu_0^{(j)})$. The difference between two such probabilities will be an eigenvalue of $\rho-\rho^\prime$, so the eigenvalues are all $\lambda_{\mathbf{n}}=p_{\mathbf{n}}-p_{\mathbf{n}}^\prime=\prod_{j=1}^M P(n_j|\mu_0^{(j)})-\prod_{j=1}^M P(n_j|\mu_0^{\prime(j)})$. The trace distance is half of the sum of the absolute values of these eigenvalues.
The question becomes: can we make anything out of these differences? The answer is marginally yes, since we specified thermal states
$$D(\rho,\rho^\prime)=\sum_{n_1,n_2,\cdots,n_M=0}^\infty\left|\prod_{j=1}^M\frac{\mu_0^{(j)n_j}}{(\mu_0^{(j)}+1)^{n_j+1}}-\prod_{j=1}^M \frac{\mu_0^{\prime(j)n_j}}{(\mu_0^{\prime(j)}+1)^{n_j+1}}\right|/2.$$ But is that helpful? The problem is that the absolute value of the difference of the products does not separate into sums over independent modes, so everything is quite convoluted.
Now let's check whether we can even simplify this for a single mode. Assuming $\mu_0^\prime<\mu_0$, we have that $p_n^\prime>p_n$ for $n>n_0$ and $p_n^\prime<p_n$ for $n<n_0$, where $n_0=\frac{\ln \frac{\mu_0+1}{\mu_0^\prime+1}}{\ln \frac{\mu_0^\prime+1}{\mu_0+1}+\ln \frac{\mu_0}{\mu_0^\prime}}$ (found from finding $n$ such that $p_n=p_n^\prime$. Rounding down as $m=\lfloor n_0\rfloor$, all four geometric series can be summed, yielding
\begin{align}
2D(\rho,\rho^\prime)=\frac{1-2\left(\frac{\mu_0^\prime}{\mu_0^\prime+1}\right)^{m+1}}{1-\frac{\mu_0^\prime}{\mu_0^\prime+1}}-\frac{1-2\left(\frac{\mu_0}{\mu_0+1}\right)^{m+1}}{1-\frac{\mu_0}{\mu_0+1}}.
\end{align}
This is highly nonlinear, as the exponent $m$ itself depends on $\mu_0$ and $\mu_0^\prime$.
To extend this to multimode states, one would need to identify regions where the probability differences are negative and positive. This would split all of the sums into small intervals where the geometric series could be computed. After adding and subtracting all of the alternating terms, one would find a final result. It is "straightforward" to extend the present calculation for a particular $p_{\mathbf{n}}$ that holds all components $n_i$ constant except for one and evaluates how the sign of $p-p^\prime$ changes; doing this when all components can change dramatically increases the number of times when the sign changes.
We can try a simplified multimode case where all of the thermal occupations are the same for a particular state: $\mu_0^{(j)}=\mu$ and $\mu_0^{\prime (j)}=\mu^\prime$. The products simplify nicely to things like $\left(\frac{\mu}{\mu+1}\right)^{N}\frac{1}{(\mu+1)^M}$ for $N=\sum_j n_j$. As before we find the crossing point $N_0=M\frac{\ln\frac{\mu+1}{\mu^\prime+1}}{\ln\frac{\mu}{\mu^\prime}+\ln\frac{\mu^\prime+1}{\mu+1}}$ and round down as $m=\lfloor N_0\rfloor$. Since all the terms agree in sign on either side of $N_0$, we can combine everything for both the cases when $\mu>\mu^\prime$ and $\mu<\mu^\prime$ by putting one absolute value outside of everything
$$2D(\rho,\rho^\prime)=\left|\sum_{N\leq m}\left(\sum_{n_1+\cdots+n_M=N}\right)\left(\left(\frac{\mu}{\mu+1}\right)^{N}\frac{1}{(\mu+1)^M}-\left(\frac{\mu^\prime}{\mu^\prime+1}\right)^{N}\frac{1}{(\mu^\prime+1)^M}\right)-\sum_{N> m}\left(\sum_{n_1+\cdots+n_M=N}\right)\left(\left(\frac{\mu}{\mu+1}\right)^{N}\frac{1}{(\mu+1)^M}-\left(\frac{\mu^\prime}{\mu^\prime+1}\right)^{N}\frac{1}{(\mu^\prime+1)^M}\right)\right|.$$ Next we simplify: how many combinations of sums of integers lead to $N$? That's like asking how may ways are there to put $N$ balls into $M$ bins, so the answer is $\binom{N+M-1}{N}$. Finally, we sum. Unfortunately, I cannot do better than converting $\sum_{N=0}^m X^N \binom{N+M-1}{N}$ into a hypergeometric function (maybe try Mathematica), so the answer will involve two differences between pairs of hypergeometric functions.
