I want to determine the vaccum expectation value of a string of creation and annihilation operators. They have a very specific form:
$$\langle \prod_{i=1}^n \hat{a}(k_i) \, \, \hat{N}_1 \prod_{j=1}^n \hat{a}^\dagger(k_j) \rangle_0,$$
which is the number operator $\hat{N}_1 = \hat{a}^\dagger(k_1) \hat{a}(k_1)$ sandwitched on the left by $n$ annhilation operators and on the right by $n$ creation operators. The subscript $0$ reminds us that it is the vacuum expectation that we are interested in. However, for any $n$ which can be arbitrarily large, there exist related expectation values that I would like to evaluate. For example, I also have
$$\langle \prod_{i=1}^n \hat{a}(k_i) \, \, \hat{a}^\dagger(k_1) \hat{N}_2 \prod_{j=2}^n \hat{a}^\dagger(k_j) \rangle_0,$$
where the number operator is positioned between $\hat{a}^\dagger(k_1)$ and $\hat{a}^\dagger(k_2)$. Note the index $j$ on the second product. This continues for all $\hat{N}_m,$ $m \in [1, n]$, i.e.
$$\langle \prod_{i=1}^n \hat{a}(k_i) \, \, \prod_{j=1}^m \hat{a}^\dagger(k_j) \hat{N}_m \prod_{l=m}^n \hat{a}^\dagger(k_l) \rangle_0.$$
Now the operators obeys the usual commutation relations, but since the arguments $k_i$ are continuous, we have Dirac delta functions:
$$\left[\hat{a}(k_i), \hat{a}^\dagger(k_j)\right] = \delta(k_i - k_j).$$
The whole expectation is in a multidimensional integral $\int dk_1 \cdots dk_n$, and so the re-expression of the expectation values in terms of Dirac delta functions will greatly simplify the problem.
My question is how can I obtain a general expression for the expectation values for arbitrary $n$? I have tried normally ordering the operator string up to $n=3$, but this becomes labourious and I haven't managed to find a general expression yet.
Can the Husimi Q-function be used to do this?
EDIT IN RESPONSE TO ANSWER BY ZEROTHEHERO
A way to solve such VEVs is to use the so-called XD-representation. However, I have been unable to calculate the above VEVs using the method. Any help with this would be appreciated!