Hubbard-Stratonovich transformation in the operator form I was reading the Chapter 14 of the textbook by Philip Phillips on Advanced Solid State Physics, when he introduced the mean-field treatment of the quantum rotor model, he used the method first used in this paper: Quantum fluctuations in two-dimensional superconductors. 
Here my question is related to the Hubbard-Stratonovich transformation (H-S transformation) used here to decouple the "spin" interaction:
$$ T_{\tau}\exp \left( \int_0^{\beta} d\tau\ \sum_k J_k \mathbf{S}_k(\tau)\cdot\mathbf{S}_k^*(\tau) \right)$$
here the $\mathbf{S}_k(\tau)$ is the operator written in the interaction picture, then he introduced the auxiliary field $\psi_k$:
$$T_{\tau}\exp \left( \int_0^{\beta} d\tau\ \sum_k J_k \mathbf{S}_k(\tau)\cdot\mathbf{S}_k^*(\tau) \right)=T_{\tau}\int D \psi_k(\tau) \ e^{-\int d\tau \sum_k \psi_k^*(\tau)\psi_k(\tau)-2\int_0^{\beta} d \tau \sum_k \psi_k(\tau)\cdot\mathbf{S}_{-k}(\tau)  } $$ 
because I learned H-S transformation in the context of path integral, where we write the partition function in terms of functional integral over the  fields ("numbers"). So here I want to ask:


*

*Mathematically if it is well defined to do the H-S transformation on top of operators?

*It seems that the imaginary time ordering operator does not affect the auxiliary field $\psi_k(\tau)$ (can be seen from later on derivations in the book), what's the reason for that? 

 A: Well-defined mathematically depends on how much one wants to be rigorous. Here, I think a physicist would be satisfied with an argument like : working with eigenvectors of ${\bf S_k}$, one shows that the formal manipulation with operators is true when working with the eigenvalues, so it is fine to do that in terms of the operators. I don't know how much effort is needed to make that satisfying for a mathematician.
Concerning the second question : the auxiliary field are (complex) numbers, so they all commute. Therefore, there isn't any ambiguity in their ordering.
By construction, in a path integral formalism, the fields are always "time-ordered". 
Update : to give a bit more details on the HS transformation for operators, I will use a toy model. Let $\hat A$ be an operator, with eigenstates $|a\rangle$ associated to the eigenvalue $a$. Then 
$$ e^{\frac12\hat A^2}|a\rangle = e^{\frac12a^2}|a\rangle,\\
=\int dx e^{-\frac12x^2+a x}|a\rangle,\\
= \int dx e^{-\frac12x^2+\hat A x}|a\rangle.$$
Thus (assuming that the eigenstates of $\hat A$ form a basis and so on), we can say that at the operator level,
$ e^{\frac12\hat A^2}=\int dx e^{-\frac12x^2+\hat A x}$ (up to a normalization constant).
