I would like to discretize the following integral operator:
$$\frac{1}{s^2}\sum_{j=1}^N\mu_j\int d\mathbf{x}d\mathbf{x}'f(\mathbf{x})f(\mathbf{x'}) \left(x_j + x'_j - 2\mu_j\right)\hat{a}^\dagger(\mathbf{x})\vert0\rangle\langle0\vert\hat{a}(\mathbf{x'}),$$
where $d\mathbf{x} = dx_1 \cdots dx_N$ and $f(\mathbf{x}) = f(x_1) \cdots f(x_N)$ is the square root of a Gaussian such that
$$f(x_i) = \frac{1}{(2\pi \sigma^2)^{1/4}} \exp\left[-\frac{(x_i - \mu_i)^2}{4 \sigma^2}\right]$$
I would like to discretize it so that it can be written in matrix form and its eigenbasis understood. I saw a similar post here, using Simpson's rule for the discretization, subsequently yielding the matrix.
Context: In quantum statistical estimation theory, the Symmetric Logarithmic derivative is an important operator. It helps to define the quantum Fisher information but more importantly for pure states, its basis defines the optimal estimation strategy which saturates the corresponding quantum Cramer-Rao Bound.
