Let $\sigma_k$ be the $k$-th elementary symmetric polynomial in $M$ variables $b_i$ and define $$\prod_{i=1}^M \frac{t b_i} {1-e^{-t b_i}} =: \sum_{m=0} t^m T_m (b)$$
Of course it's easy to express $T_m(b)$ in terms of $b$ and Bernoulli, but is there a closed formula to write $T_m(b)$ as $T_m(\sigma_1,\ldots,\sigma_m)$, in terms of Bernoulli and $\sigma_1,\ldots,\sigma_m$, where $\sigma_k = \sigma_k(b)$?
For example, copying from Wikipedia, the first few terms should be $$ 1 + \frac{\sigma_1}2 + \frac{\sigma_1^2+\sigma_2}{12} + \frac{\sigma_1\sigma_2}{24} + \frac{-\sigma_1^4 + 4\sigma_1^2\sigma_2 + \sigma_1\sigma_3 + 3\sigma_2^2 − \sigma_4}{720} + \ldots$$