Normal ordered exponential of one-body operators 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.
 A: Hints: First try to show it for a single bosonic oscillator (for fermions this was done by the OP already). To this end, define the following functions:
\begin{align}
f(M)&:=\exp{a^\dagger a M} \tag{1} \\
L(M)&:=N[\exp{a^\dagger a (e^{M}-1)}] \quad  \tag{2}.
\end{align}
Then show $f(0)=L(0)=\mathbb I$ and that $f$ and $L$ satisfy the same differential equation, which in turn implies $f(M)=L(M)$, proving the claim.
The generalizations for $N>1$ and fermions are left to you.
Here are some useful relations you can use/find/prove:
\begin{align}
e^{-a^\dagger a M}\, a\, e^{a^\dagger aM} &=  e^M\, a \tag{3}\\
N[(a^\dagger a)^n] &= (a^\dagger)^n a^n \tag{4}\\
f^\prime(M) &= a^\dagger a\, f(M) \overset{(3)}{=} e^M a^\dagger\, f(M)\, a \tag{5}\quad .
\end{align}
A: Let us prove OP's claim for a single bosonic mode:

Proposition:
$$ e^{ta^{\dagger}a}~=~:e^{(e^t-1)a^{\dagger}a}: \qquad t~\in~\mathbb{C}.\tag{A}$$

Sketched proof of eq. (A): Let's call the LHS for $U(t)$ and the RHS for $V(t)$. Both sides can be written as a function of the operator $n=a^{\dagger}a$ without the use of $a$ and $a^{\dagger}$. They satisfy the same first-order ODE:
$$U^{\prime}(t)~=~a^{\dagger}a e^{tn}~=~a^{\dagger} e^{t(n+1)}a~=~e^t a^{\dagger} U(t)a, \tag{B}$$
$$V^{\prime}(t)~=~e^t a^{\dagger} V(t)a, \tag{C}$$
with the same initial condition $U(0)={\bf 1}=V(0)$. Hence they must be equal. $\Box$
