Let $\{a_i\}_{i=1}^N$ be a set of annihilation operators (they are either all bosons, or all fermions) satisfying the canonical commutation or anti-commutation relation. In the book Quantum Theory of Finite Systems by Blaizot and Ripka, Problem 1.6 claims that (summation over repeated indices is implied)
$$ \exp(a^\dagger_i M_{ij} a_j) = N[\exp(a^\dagger_i (e^M-1)_{ij} a_j)] \tag{1} $$
where $M_{ij}$ is an $N \times N$ complex matrix, and $N[A]$ puts creation operators in $A$ to the left, treating all $a^\dagger, a$ in the argument as commuting or anti-commuting numbers. For example, with $\eta = +1$ for bosons, and $-1$ for fermions, we have
$$ N[a_4 a^\dagger_2 a_1 a^\dagger_3] = \eta^{1 + 2} a^\dagger_2 a^\dagger_3 a_4 a_1 = \eta a^\dagger_2 a^\dagger_3 a_4 a_1 $$
I tried to prove eq. (1) by series expansion and comparing terms, but the expansion soon becomes rather complicated. I would appreciate it if someone can provide an elegant and clean proof.
My current attempt: For fermions the exponential function can be greatly simplified. Below I give a proof for fermions when $N = 1$, so that $M$ reduced to a complex number.
The RHS (right hand side) of eq. (1) now actually means
$$ \begin{align*} \text{RHS} &= N[\exp[(e^{M}-1) a^\dagger a]] \\ &= 1 + \sum_{n=1}^\infty \frac{(e^{M}-1)^n}{n!} N\left[(a^\dagger a)^n\right] \end{align*} $$
Normal ordering gives:
$$ \begin{align*} N\left[(a^\dagger a)^n\right] &= N[a^\dagger a a^\dagger a \cdots a^\dagger a] \\ &= \eta^{1 + \cdots + (n-1)} a^{\dagger n} a^n \\ &= \eta^{n(n-1)/2} a^{\dagger n} a^n \end{align*} $$
For fermions, $a^n = 0$ for $n \ge 2$, which is the key to simplify the exponential function:
$$ \begin{align*} \text{RHS} &= 1 + (e^M - 1) a^\dagger a \end{align*} $$
Meanwhile,
$$ \begin{align*} \text{LHS} &= \exp(M a^\dagger a) = 1 + \sum_{n=1}^\infty \frac{M^n}{n!} (a^\dagger a)^n \end{align*} $$
But with $a a^\dagger = 1 - a^\dagger a$, we notice that
$$ \begin{align*} (a^\dagger a)^2 &= a^\dagger a a^\dagger a = a^\dagger (1 - a^\dagger a) a \\ &= a^\dagger a - \underbrace{ a^{\dagger 2} a^2 }_{= 0} = a^\dagger a \end{align*} $$
which further leads to $(a^\dagger a)^n = a^\dagger a$ for any $n \ge 1$. Therefore
$$ \begin{align*} \text{LHS} &= 1 + \bigg[ \sum_{n=1}^\infty \frac{M^n}{n!} \bigg] a^\dagger a \\ &= 1 + (e^M - 1) a^\dagger a = \text{RHS} \end{align*} $$
But obviously things will be complicated for bosons, since the $a$ operator is no longer nilpotent.