Gaussian Grassmann integral with complex/bosonic source term I'm interested in solving the following multi-dimensional integral
$$
\int d \theta d \bar{\theta} e^{-\bar{\theta}M \theta +\Lambda \theta + \bar{\theta} J }
$$
where $\theta$ is a $N$-dimensional vector of odd Grassmann numbers while $M,\Lambda,J$ are composed of even Grassmann numbers (different from $\theta$). I have looked everywhere but I can only find situations where $\Lambda,J$ are also odd Grassmann variables. In that case, we can do a simple shift of variables such that the new variables also have a well-defined odd-parity. This is not the case for even $\Lambda,J$.
A simpler case of this integral is that of a single Grassmann number and $M,\Lambda,J$ being complex matrices. I've tried to solve by expanding the exponential and ignoring the higher-order terms $\mathcal{O}(\bar{\theta}\theta\bar{\theta}\theta)$
$$
\int d \theta d \bar{\theta} e^{-\bar{\theta}  M \theta +\Lambda \theta + \bar{\theta} J }=\int d \theta d \bar{\theta} (1-\bar{\theta } M \theta +\Lambda \theta + \bar{\theta} J +\frac{1}{2} (\Lambda \theta \bar{\theta} J + \bar{\theta} J \Lambda \theta))=\int d \theta d \bar{\theta} (-M \bar{\theta }\theta +\frac{\Lambda J}{2} (\theta \bar{\theta} + \bar{\theta} \theta)).
$$
but since $\theta$ are even numbers we have $(\theta \bar{\theta} + \bar{\theta} \theta)=0$ and the final result does not depend on $J,\Lambda$ which should be wrong. Could you help me understand what is wrong? Thank you.
 A: OP's manipulations are correct for $N=1$ if we define the exponential function via its Taylor series, cf. Ref. 1. Note that the argument of the exponential function is usually Grassmann-even in physical theories.
References:

*

*Bryce DeWitt, Supermanifolds, Cambridge Univ. Press, 1992; eq. (1.1.6).

A: I believe I have constructed a correct answer that is independent of the parity of the shift variable but the result does not agree with the 1-dimensional case discussed in the question.
$$ \int d\theta d\bar{\theta}e^{-\bar{\theta}M\theta+\Lambda\bar{\theta}+\theta J} \\=\int d\theta d\bar{\theta}e^{-\tilde{\bar{\theta}}M\tilde{\theta}+\Lambda M^{-1}J}
 \\=\int d\theta d\bar{\theta}\sum_{n}\frac{\left(\Lambda M^{-1}J-\tilde{\bar{\theta}}M\tilde{\theta}\right)^{n}}{n!}
 \\=\int d\theta d\bar{\theta}\sum_{n}\frac{\text{Permutations}\left\{ \Lambda M^{-1}J,-\tilde{\bar{\theta}}M\tilde{\theta}\right\} }{n!}$$
Where we redefined the variables
$$\tilde{\theta}=\theta-M^{-1}J\qquad\text{and}\qquad\tilde{\bar{\theta}}=\bar{\theta}-\Lambda M^{-1}.$$
By Permutations$\{A,B\}$ I mean all the terms $AAA...,BAA...,ABA...,...BBB$ without exchanging any terms since $\tilde{\theta}$ does not have a well defined parity and we want to avoid it. Next we recall that only the highest order terms in $\theta,\bar{\theta}$ can remain and those correspond to replacing $\tilde{\theta}\rightarrow \theta$. But since $\Lambda M^{-1}J$ and $\bar{\theta}M\theta$ both have even parity, we can now apply binomial theorem to get:
$$ \int d\theta d\bar{\theta}\sum_{n}\frac{\text{Permutations}\left\{ \Lambda M^{-1}J,-\bar{\theta}M\theta\right\} }{n!}+\text{lower order terms in $\theta$}
 \\=\int d\theta d\bar{\theta}\sum_{n=0}\sum_{k=0}^{n}\left(\begin{array}{c}
n\\
k
\end{array}\right)\frac{\left(\Lambda M^{-1}J\right)^{n-k}\left(-\bar{\theta}M\theta\right)^{k}}{n!}
 \\=\int d\theta d\bar{\theta}\sum_{n=d}\left(\begin{array}{c}
n\\
d
\end{array}\right)\frac{\left(\Lambda M^{-1}J\right)^{n-d}\left(-\bar{\theta}M\theta\right)^{d}}{n!}
 \\=\frac{e^{\Lambda M^{-1}J}}{d!}\int d\theta d\bar{\theta}\left(-\bar{\theta}M\theta\right)^{d}
 \\=e^{\Lambda M^{-1}J}\det M$$
