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.
Best Answer
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}