Correlated Gaussian method in Quantum Mechanics

Explicitly correlated Gaussian functions have been extensively used in quantum-mechanical variational calculations in atomic, molecular, and nuclear physics. This book is an attempt to collect the relevant information about this established tool in computational quantum mechanics.

Introduction

Often we describe bound-state and scattering problems in nuclear and atomic physics through the Shrödinger equation. Unfortunately modern quantum physics offers problems that we can not solve analytically. Luckily the availability of powerful computers is shifting the emphasis from the analytical computation of the solution toward numerical analysis. During the last century numerous methods were developed in order to approximate solutions numerically, e.g. Monte-Carlo simulation, Hypershperical expansion, variational methods with different trial wave functions etc.

In this section we discuss variational method with trial function in the form of correlated Gaussians which is widely used in the modern physics. Mathematically it is based on the Ritz theorem, that states that for an arbitrary function Ψ from the state space the expectation value of the Hamiltonian (<Ψ|H|Ψ>/<Ψ|Ψ>) is larger then the ground state energy. So choosing different trial wave functions and calculating mean values of the Hamiltonian for this functions allows us to get an upper bound for the ground state energy.

Example

To show the idea of the method we consider two particles in 1 dimension interacting through the oscillator potential

${\bigg [}-{\frac {1}{2}}{\frac {\partial ^{2}}{\partial {x}_{1}}}-{\frac {1}{2}}{\frac {\partial ^{2}}{\partial {x}_{2}}}+{\frac {1}{2}}(x_{1}-x_{2})^{2}{\bigg ]}\phi =E\phi \;.$

It is really simple textbook problem with the ground state solution

$\phi _{0}=({\frac {\sqrt {2}}{\pi }})^{1/4}\exp(-{\frac {1}{2{\sqrt {2}}}}(x_{1}-x_{2})^{2})\;,E_{0}={\frac {1}{\sqrt {2}}}\;,$

where we assumed for simplicity, that the total momentum is equal to zero.

To show how the method works we choose the trial wave function in the form of just one gaussian

$f_{tr}=({\frac {4\alpha }{\pi }})^{1/4}\exp(-\alpha (x_{1}-x_{2})^{2})\;,$

where we have only one real positive parameter which has to minimize the energy. The idea of the method is to pick this parameter stochastically using just generator of the real numbers. We find out that independently of the seed after 50 attempts we found value of α that gives the ground state energy with 5 significant digits.

Convergence of the expectation value of the Hamiltonian (different curves corresponds to different starting seeds for the generator of the random numbers).

Of course it's really simple example and we can establish energy with high precision just because we are working in the space that contains ground state wave function of the Hamiltonian.

In order to establish excited states it's not enough to use just one Gaussian so we pick trial wave function in more general form

$f_{tr}=\sum _{i=1}^{N}c_{i}\exp(-\alpha _{i}(x_{1}-x_{2}-s_{i})^{2})\;.$

As above we pick parameters $\alpha _{i},s_{i}$ stochastically and then determine linear parameters c_i demanding minimal expectation value of the Hamiltonian. Using N=25 and one random set of the parameters $\alpha _{i},s_{i}$ (we assume that $0<\alpha _{i}<5,-1<s_{i}<1$) we get first 6 eigenstates with 5 significant digits.

From this simple example we learn that we are able to approximate solution of the Schrödinger equation without any preliminary knowledge about the system using only random search. The main problem is to estimate how good this approximation is.

The method

Hamiltonian under consideration

We shall consider a non-relativistic quantum mechanical $N$-particle (N-body) system described by a Hamiltonian

${\hat {H}}=\sum _{i=1}^{N}\left(-{\frac {\hbar ^{2}}{2m_{i}}}\right){\frac {\partial ^{2}}{\partial {\vec {r}}_{i}}}+\sum _{i=1}^{N}V_{i}({\vec {r}}_{i})+\sum _{i<j}^{N}V_{ij}({\vec {r}}_{i}-{\vec {r}}_{j})\;,$

where $m_{i}$ and ${\vec {r}}_{i}$ are the mass and the coordinate of particle number $i$; the first term is the kinetic energy operator; the second term is the one-body force, like the external field (often an oscillator, $V_{i}({\vec {r}}_{i})={\frac {m_{i}\omega ^{2}}{2}}{\vec {r}}_{i}^{2}$), for example Magneto-optical trap; and the third term is the two-body forces, the inter-particle interactions.

This Hamiltonian can describe a system of atoms in a trap or a nucleus.

In general one can also consider spin-dependent, three-body, non-local and other types of interactions which we shall introduce when needed.

Basis expansion of the Schrödinger equation

We are going to solve the Schrödinger equation

${\hat {H}}|\psi \rangle =\epsilon |\psi \rangle \;,$

where ${\hat {H}}$ is the Hamiltonian of a quantum few-body system, $|\psi \rangle$ and $\epsilon$ are the eigenfunction and the eigenvalue to be found.

We shall expand the wave-function $|\psi \rangle$ in terms of a set of basis functions $|i\rangle _{i=1\dots n}$,

$|\psi \rangle =\sum _{i=1}^{n}c_{i}|i\rangle \;.$

Inserting the expansion into the Schrödinger equation and multiplying from the left with $\langle k|_{1\leq k\leq n}$ gives

$\sum _{i=1}^{n}\langle k|{\hat {H}}|i\rangle c_{i}=\epsilon \sum _{i=1}^{n}\langle k|i\rangle c_{i}\;.$

Or, in the matrix notation

${\cal {H}}c=\epsilon {\cal {N}}c\;,$

where ${\cal {H}}$ and ${\cal {N}}$ are correspondingly the Hamiltonian and the overlap matrices with the matrix elements

${\cal {H}}_{ki}=\langle k|{\hat {H}}|i\rangle \,,\;{\cal {N}}_{ki}=\langle k|i\rangle \;.$

The matrix equation ${\cal {H}}c=\epsilon {\cal {N}}c$ is called the generalized eigenvalue problem. There exist established routines to solve this problem, for example the function "eig" in Octave, or the function "gsl_eigen_gensymm" in the GNU Scientific Library.

If the basis function are orthogonal and normalized the overlap matrix equals the unity matrix, ${\cal {N}}_{ij}=\delta _{ij}$, and the generalized eigenvalue problem reduces to the ordinary eigenvalue problem, ${\cal {H}}c=\epsilon c$.

We shall use Gaussians --- which are not orthogonal --- as basis functions, therefore we shall deal with the generalized eigenvalue problem.

Gaussians as basis functions

We shall use the so-called Correlated Gaussians (or Explicitly Correlated Gaussians) as the basis functions. For a system of $N$particles with coordinates ${\vec {r}}_{i}\,,\;i=1\dots N$, the Correlated Gaussian is defined as

$g({\vec {r}}_{1},\dots ,{\vec {r}}_{N})\doteq \exp \left(-\sum _{i,j=1}^{N}A_{ij}{\vec {r}}_{i}\cdot {\vec {r}}_{j}+\sum _{i=1}^{N}{\vec {s}}_{i}\cdot {\vec {r}}_{i}\right)$,

where ${\vec {r}}_{i}\cdot {\vec {r}}_{j}$ denotes the dot-product of the two vectors; and where $A$, a symmetric positive-defined matrix, and ${\vec {s}}_{i}|_{i=1,\dots ,N}$, the shift-vectors, are (cleverly chosen) parameters of the Gaussian.

In matrix notation,

$g({\bf {r}})\doteq \exp \left(-{\bf {r}}^{\mathsf {T}}A{\bf {r}}+\mathbf {s} ^{\mathsf {T}}\mathbf {r} \right)$,

where ${\bf {r}}$ is the column of the coordinates ${\vec {r}}_{i}$ and $\mathbf {s}$ is the column of the shift-vectors ${\vec {s}}_{i}$,

${\bf {r}}\doteq \left({\begin{matrix}{\vec {r}}_{1}\\\vdots \\{\vec {r}}_{N}\end{matrix}}\right),\;{\bf {s}}\doteq \left({\begin{matrix}{\vec {s}}_{1}\\\vdots \\{\vec {s}}_{N}\end{matrix}}\right)$,

and

${\bf {r}}^{\mathsf {T}}A{\bf {r}}+\mathbf {s} ^{\mathsf {T}}\mathbf {r} \doteq \sum _{ij}{\vec {r}}_{i}\cdot A_{ij}{\vec {r}}_{j}+\sum _{i}{\vec {s}}_{i}\cdot {\vec {r}}_{i}$.

Matrix elements with Gaussians

We represent the few-body wave-function $\Psi ({\vec {r}}_{1},\dots ,{\vec {r}}_{N})$as a linear combination of $n$ correlated Gaussians,

$\Psi ({\vec {r}}_{1},\dots ,{\vec {r}}_{N})=\sum _{i=1}^{n}c_{i}g_{i}({\vec {r}}_{1},\dots ,{\vec {r}}_{N})\;.$

We shall discuss the choice of matrices $A$ in these Gaussians later and here only calculate the matrix elements.

Overlap

The overlap,

$\langle g'|g\rangle \doteq \int d^{3}{\vec {r}}_{1}\dots d^{3}{\vec {r}}_{N}\exp \left(-{\bf {r}}^{\mathsf {T}}(A'+A){\bf {r}}+(\mathbf {s} '+\mathbf {s} )^{\mathsf {T}}\mathbf {r} \right)\equiv \int d\mathbf {r} \exp \left(-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} \right)\;,$

(where $d\mathbf {r} \doteq d^{3}{\vec {r}}_{1}\dots d^{3}{\vec {r}}_{N}$, $B\doteq A'+A$ and $\mathbf {v} \doteq \mathbf {s} '+\mathbf {s}$) can be calculated making an orthogonal transformation, ${\bf {r}}\rightarrow {\bf {x}}=Q^{\mathsf {T}}{\bf {r}}$, where $Q^{\mathsf {T}}Q=1$, to the basis where the matrix $B$ is diagonal,

${\begin{aligned}\langle g'|g\rangle &=\int d\mathbf {x} \exp \left(-{\bf {x}}^{\mathsf {T}}B{\bf {x}}+\mathbf {v} ^{\mathsf {T}}\mathbf {x} \right)=\int d\mathbf {x} \exp \left(-\sum _{i=1}^{N}{\vec {x}}_{i}\cdot B_{ii}{\vec {x}}_{i}+\sum _{i=1}^{N}{\vec {v}}_{i}\cdot {\vec {x}}_{i}\right)&\\&=\prod _{i=1}^{N}\int d^{3}{\vec {x}}_{i}\exp \left(-{\vec {x}}_{i}\cdot B_{ii}{\vec {x}}_{i}+{\vec {v}}_{i}\cdot {\vec {x}}_{i}\right)=\prod _{i=1}^{N}\exp \left({\frac {1}{4B_{ii}}}{\vec {v}}_{i}^{2}\right)\left({\frac {\pi }{B_{ii}}}\right)^{3/2}\\&=e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\;.&\end{aligned}}$

Finally,

$\langle g'|g\rangle =e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\,,\;\mathrm {where} \;B\doteq A'+A\;\mathrm {and} \;\mathbf {v} \doteq \mathbf {s} '+\mathbf {s} \;.$

Kinetic energy

Let us consider a more general form of the kinetic energy operator,

${\hat {K}}=-\sum _{i,j=1}^{N}{\frac {\partial }{\partial {\vec {r}}_{i}}}\Lambda _{ij}{\frac {\partial }{\partial {\vec {r}}_{j}}}\equiv -{\frac {\partial }{\partial {\bf {r}}}}\Lambda {\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}\;,$

where $\Lambda$ is a constant matrix, for example, for the Hamiltonian above $\Lambda _{ij}={\frac {\hbar ^{2}}{2m_{i}}}\delta _{ij}\;.$

The matrix element is given as

$\langle g'|{\hat {K}}|g\rangle =\int d\mathbf {r} e^{-{\bf {r}}^{\mathsf {T}}A'{\bf {r}}+\mathbf {s} '^{\mathsf {T}}\mathbf {r} }\left(-\sum _{i,j=1}^{N}{\frac {\partial }{\partial {\vec {r}}_{i}}}\Lambda _{ij}{\frac {\partial }{\partial {\vec {r}}_{j}}}\right)e^{-{\bf {r}}^{\mathsf {T}}A{\bf {r}}+\mathbf {s} ^{\mathsf {T}}\mathbf {r} }\;.$

Integrating by parts with respect to the left derivative gives

${\begin{aligned}\langle g'|{\hat {K}}|g\rangle &=\int d\mathbf {r} \left({\frac {\partial }{\partial \mathbf {r} }}e^{-{\bf {r}}^{\mathsf {T}}A'{\bf {r}}+\mathbf {s} '^{\mathsf {T}}\mathbf {r} }\right)\Lambda \left({\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}e^{-{\bf {r}}^{\mathsf {T}}A{\bf {r}}+\mathbf {s} ^{\mathsf {T}}\mathbf {r} }\right)\\&=\int d\mathbf {r} e^{-{\bf {r}}^{\mathsf {T}}B{\bf {r}}+\mathbf {v} ^{\mathsf {T}}\mathbf {r} }(-2\mathbf {r} ^{\mathsf {T}}A'+\mathbf {s} '^{\mathsf {T}})\Lambda (-2A\mathbf {r} +\mathbf {s} )\\&=\int d\mathbf {r} e^{-{\bf {r}}^{\mathsf {T}}B{\bf {r}}+\mathbf {v} ^{\mathsf {T}}\mathbf {r} }(4\mathbf {r} ^{\mathsf {T}}A'\Lambda A\mathbf {r} -2\mathbf {r} ^{\mathsf {T}}A'\Lambda \mathbf {s} -2\mathbf {s} '^{\mathsf {T}}\Lambda A\mathbf {r} +\mathbf {s} '^{\mathsf {T}}\Lambda \mathbf {s} )\;.\end{aligned}}$

Calculating term after term,

${\begin{aligned}\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} }\mathbf {r} ^{\mathsf {T}}F\mathbf {r} &={\frac {\partial }{\partial \mathbf {v} }}F{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} }=\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\left({\frac {\partial }{\partial \mathbf {v} }}F{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\\&=\left({\frac {3}{2}}\;\mathrm {trace} (FB^{-1})+\mathbf {u} ^{\mathsf {T}}F\mathbf {u} \right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\;,\end{aligned}}$

${\begin{aligned}\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} }\mathbf {a} ^{\mathsf {T}}\mathbf {r} &=\mathbf {a} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} }=\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\left(\mathbf {a} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\\&=\left(\mathbf {a} ^{\mathsf {T}}\mathbf {u} \right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\;,\end{aligned}}$

where $\mathbf {u} \doteq {\frac {1}{2}}B^{-1}\mathbf {v}$.

Finally,

$\left\langle g'\left|-{\frac {\partial }{\partial {\bf {r}}}}\Lambda {\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}\right|g\right\rangle =\left(6\,\mathrm {trace} (A'\Lambda AB^{-1})+(2A'\mathbf {u} -\mathbf {s} ')^{\mathsf {T}}\Lambda (2A\mathbf {u} -\mathbf {s} )\right)e^{{\frac {1}{2}}\mathbf {v} ^{\mathsf {T}}\mathbf {u} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\,,\;B\doteq A'+A,\;\mathbf {u} \doteq {\frac {1}{2}}B^{-1}\mathbf {v} .$

Potential energy

The matrix element of the interaction potential between, say, particles $i$ and $j$,

$\langle g'|V({\vec {r}}_{i}-{\vec {r}}_{j})|g\rangle$,

can be written in a more general form,

$\langle g'|V(w^{\mathsf {T}}\mathbf {r} )|g\rangle$,

where $w$ is a size-N vector with all components equal zero except for $w_{i}=1$ and $w_{j}=-1$.

A one-body interaction, $V({\vec {r}}_{i})$, has the same form, $V(w^{T}{\bf {r}})$, where $w_{i}=1$ and all other components equal zero.

Central potential

For a central potential in a Gaussian form,

$V({\vec {r}}_{i})\doteq Se^{-\alpha {\vec {r}}_{i}^{2}}=Se^{-\alpha \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} },\;\mathrm {where} \;w_{i}=1,\;w_{k\neq i}=0\;,$

the matrix element between shifted Gaussians can be calculated straightforwardly,

$\langle g'|e^{-\alpha \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }|g\rangle =\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}(B+\alpha ww^{\mathsf {T}})\mathbf {r} +\mathbf {v} ^{\mathsf {T}}\mathbf {r} }=e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}$

where $B'=B+\alpha ww^{\mathsf {T}}$.

The rank-1 updates of the determinant, $\det(B+\alpha ww^{\mathsf {T}})$, and the matrix inverse, $(B+\alpha ww^{\mathsf {T}})^{-1}$, can be efficiently calculated using the following formulae,

$\det(B+uv^{\mathsf {T}})=(1+v^{\mathsf {T}}B^{-1}u)\,\det(B)\,,$

and (Sherman-Morrison)

$(B+uv^{T})^{-1}=B^{-1}-{B^{-1}uv^{T}B^{-1} \over 1+v^{T}B^{-1}u}.$

For a general form-factor central potential one way to calculate the matrix element is through the Fourier transform of the potential, $V({\vec {r}})=\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}{\mathfrak {f}}({\vec {k}})e^{i{\vec {k}}\cdot {\vec {r}}}$,

$\langle g'|V(w^{\mathsf {T}}\mathbf {r} )|g\rangle =\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}{\mathfrak {f}}({\vec {k}})\int d{\bf {r}}e^{-{\bf {r}}^{\mathsf {T}}B{\bf {r}}+(\mathbf {v} +i{\vec {k}}w)^{\mathsf {T}}{\bf {r}}}=\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}{\mathfrak {f}}({\vec {k}})\;e^{{\frac {1}{4}}(\mathbf {v} +i{\vec {k}}w)^{\mathsf {T}}B^{-1}(\mathbf {v} +i{\vec {k}}w)}\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}$.

Thus

$\langle g'|V(w^{\mathsf {T}}\mathbf {r} )|g\rangle =e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}{\mathfrak {f}}({\vec {k}})e^{-\alpha {\vec {k}}^{2}+i{\vec {k}}{\vec {q}}}\,,$

where ${\mathfrak {f}}({\vec {k}})$ is the Fourier transform of the potential $V({\vec {r}})$, $\alpha \doteq {\frac {1}{4}}w^{T}B^{-1}w$ and ${\vec {q}}\doteq {\frac {1}{2}}w^{\mathsf {T}}B^{-1}\mathbf {v}$.

The last integral can be also written via the potential itself,

${\begin{aligned}\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}{\mathfrak {f}}({\vec {k}})e^{-\alpha {\vec {k}}^{2}+i{\vec {k}}{\vec {q}}}&=\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}\int d^{3}{\vec {r}}V({\vec {r}})e^{-\alpha {\vec {k}}^{2}-i{\vec {k}}({\vec {r}}-{\vec {q}})}=\int d^{3}{\vec {r}}V({\vec {r}})e^{-{\frac {({\vec {r}}-{\vec {q}})^{2}}{4\alpha }}}\int {\frac {d^{3}{\vec {k}}}{(2\pi )^{3}}}e^{-\alpha \left({\vec {k}}+i{\frac {{\vec {r}}-{\vec {q}}}{2\alpha }}\right)^{2}}=\int d^{3}{\vec {r}}V({\vec {r}})e^{-{\frac {({\vec {r}}-{\vec {q}})^{2}}{4\alpha }}}{\frac {\pi ^{3/2}}{(2\pi )^{3}\alpha ^{3/2}}}\\&=\left({\frac {1}{4\alpha \pi }}\right)^{3/2}\int d^{3}{\vec {r}}V({\vec {r}})e^{-{\frac {({\vec {r}}-{\vec {q}})^{2}}{4\alpha }}}=\left({\frac {\beta }{\pi }}\right)^{3/2}\int d^{3}{\vec {r}}V({\vec {r}})e^{-\beta ({\vec {r}}-{\vec {q}})^{2}}\;,\end{aligned}}$

where $\beta \doteq {\frac {1}{4\alpha }}=(w^{T}B^{-1}w)^{-1}$.

Finally,

$\langle g'|V(w^{\mathsf {T}}\mathbf {r} )|g\rangle =e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B)}}\right)^{3/2}\left({\frac {\beta }{\pi }}\right)^{3/2}\int d^{3}{\vec {r}}V({\vec {r}})e^{-\beta ({\vec {r}}-{\vec {q}})^{2}}\,,$

where $\beta =(w^{T}B^{-1}w)^{-1}$ and ${\vec {q}}\doteq {\frac {1}{2}}w^{\mathsf {T}}B^{-1}\mathbf {v}$.

Here are integrals for some popular potentials,

• Gaussian,$\int d^{3}{\vec {r}}\;e^{-\gamma {\vec {r}}^{2}}\;e^{-\beta ({\vec {r}}-{\vec {q}})^{2}}=e^{-{\frac {\gamma \beta }{\gamma +\beta }}{\vec {q}}^{2}}\left({\frac {\pi }{\gamma +\beta }}\right)^{3/2}$;

• Coulomb, $\int d^{3}{\vec {r}}\;{\frac {1}{r}}\;e^{-\beta ({\vec {r}}-{\vec {q}})^{2}}=2\pi \int _{0}^{\infty }rdr\int _{-1}^{1}d\cos \theta \;e^{-\beta r^{2}+2\beta rq\cos \theta -\beta q^{2}}={\frac {2\pi }{\beta }}{\frac {1}{2q}}\int _{0}^{\infty }dr\left(e^{-\beta (r-q)^{2}}-e^{-\beta (r+q)^{2}}\right)={\frac {2\pi }{\beta }}{\frac {\sqrt {\pi }}{2}}{\frac {\mathrm {erf} ({\sqrt {\beta }}q)}{{\sqrt {\beta }}q}}$;

• Oscillator, $\int d^{3}{\vec {r}}\;{\vec {r}}^{2}\;e^{-\beta ({\vec {r}}-{\vec {q}})^{2}}=\beta ^{-5/2}\int d^{3}{\vec {x}}\left({\vec {x}}+{\sqrt {\beta }}{\vec {q}}\right)^{2}e^{-{\vec {x}}^{2}}={\frac {\pi ^{3/2}}{\beta ^{5/2}}}\left({\frac {3}{2}}+\beta {\vec {q}}^{2}\right)$.

Tensor potential

In nuclear physics the tensor potential between two nucleons has the form

${\hat {V}}_{\mathsf {tensor}}=f(r)({\vec {\sigma }}_{1}\cdot {\vec {r}})({\vec {\sigma }}_{2}\cdot {\vec {r}})$

where $f(r)$ is the form-factor; ${\vec {r}}={\vec {r}}_{1}-{\vec {r}}_{2}$; ${\vec {r}}_{1}$ and ${\vec {r}}_{2}$ are the coordinates of the nucleons; $\sigma _{1}$, $\sigma _{2}$ are the Pauli matrices related to the spins, ${\vec {s}}_{1}$ and ${\vec {s}}_{1}$, of the two nucleons: ${\vec {s}}_{1}={\frac {1}{2}}{\vec {\sigma }}_{1}$, ${\vec {s}}_{2}={\frac {1}{2}}{\vec {\sigma }}_{2}$.

One often adds the term $-f(r){\frac {1}{3}}\sigma _{1}\cdot \sigma _{2}$ to make sure that the potential has no central component (that is, the average of the potential over all directions is zero). Without this extra term the above tensor potential has a central spin-spin component ${\frac {1}{3}}f(r)(\sigma _{1}\cdot \sigma _{2}$).

Again introducing the column $w$ with $w_{1}=1$, $w_{2}=-1$, $w_{i\neq 1,2}=0$, and the vector-columns $\mathbf {a} ={\vec {\sigma }}_{1}w$, $\mathbf {b} ={\vec {\sigma }}_{2}w$, the potential can be written in a convenient general form,

${\hat {V}}_{\mathsf {tensor}}=f(w^{\mathsf {T}}\mathbf {r} )(\mathbf {a} ^{\mathsf {T}}\mathbf {r} )(\mathbf {b} ^{\mathsf {T}}\mathbf {r} )$

The matrix element of this operator between shifted Gaussians is given as

$\langle g'|{\hat {V}}_{\mathsf {tensor}}|g\rangle =\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} -\mathbf {v} ^{\mathsf {T}}\mathbf {r} }f(w^{\mathsf {T}}\mathbf {r} )(\mathbf {a} ^{\mathsf {T}}\mathbf {r} )(\mathbf {b} ^{\mathsf {T}}\mathbf {r} )=\left(\mathbf {a} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)\left(\mathbf {b} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} -\mathbf {v} ^{\mathsf {T}}\mathbf {r} }f(w^{\mathsf {T}}\mathbf {r} )\;.$

This can be calculated analytically for a Gaussian form-factor, $f(w^{\mathsf {T}}\mathbf {r} )=e^{-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }$,

${\begin{aligned}\left(\mathbf {a} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)\left(\mathbf {b} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)\int d\mathbf {r} e^{-\mathbf {r} ^{\mathsf {T}}B\mathbf {r} -\mathbf {v} ^{\mathsf {T}}\mathbf {r} }e^{-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }&=\left(\mathbf {a} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)\left(\mathbf {b} ^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {v} ^{\mathsf {T}}}}\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}\\&=\left({\frac {1}{2}}\mathbf {a} ^{\mathsf {T}}B'^{-1}\mathbf {b} +{\frac {1}{4}}(\mathbf {a} ^{\mathsf {T}}B'^{-1}\mathbf {v} )(\mathbf {b} ^{\mathsf {T}}B'^{-1}\mathbf {v} )\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}\;,\end{aligned}}$

where $B'=B+\gamma ww^{\mathsf {T}}$.

Finally, for a Gaussian tensor potential,

$\left\langle g'\left|e^{-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }(\mathbf {a} ^{\mathsf {T}}\mathbf {r} )(\mathbf {b} ^{\mathsf {T}}\mathbf {r} )\right|g\right\rangle =\left({\frac {1}{2}}\mathbf {a} ^{\mathsf {T}}B'^{-1}\mathbf {b} +{\frac {1}{4}}(\mathbf {a} ^{\mathsf {T}}B'^{-1}\mathbf {v} )(\mathbf {b} ^{\mathsf {T}}B'^{-1}\mathbf {v} )\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}\;,$

where $B'=B+\gamma ww^{\mathsf {T}}$.

Spin-orbit potential

The spin-orbit potential between two nucleons --- with coordinates ${\vec {r}}_{1}$ and ${\vec {r}}_{2}$ and spins ${\frac {1}{2}}{\vec {\sigma }}_{1}$ and ${\frac {1}{2}}{\vec {\sigma }}_{2}$ --- is usually written in the form

$V_{\mathsf {spin-orbit}}=f(r)\left({\vec {S}}\cdot {\vec {L}}\right)$

where ${\vec {r}}={\vec {r}}_{1}-{\vec {r}}_{2}$; ${\vec {S}}$ is the total spin of the two nucleons,

${\vec {S}}={\frac {1}{2}}({\vec {\sigma }}_{1}+{\vec {\sigma }}_{2})\,,$

and ${\vec {L}}$ is the relative orbital momentum between the two nucleons

${\vec {L}}=({\vec {r}}_{1}-{\vec {r}}_{2})\times {\frac {-i}{2}}\left({\frac {\partial }{\partial {\vec {r}}_{1}}}-{\frac {\partial }{\partial {\vec {r}}_{2}}}\right)\,.$

where $\times$ denotes vector-product of two vectors.

The orbital momentum can be rewritten --- using the column $w$ with $w_{1}=1$, $w_{2}=-1$, $w_{i\neq 1,2}=0$ --- in a convenient general form,

${\vec {L}}={\frac {-i}{2}}\left(w^{\mathsf {T}}\mathbf {r} \times w^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}\right)\,.$

For a Gaussian form-factor, $f(w^{\mathsf {T}}\mathbf {r} )=\exp(-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} )$, the corresponding matrix element between shifted Gaussians can be calculated analytically,

${\begin{aligned}\left\langle g'\left|e^{-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }\left(w^{\mathsf {T}}\mathbf {r} \times w^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}\right)\right|g\right\rangle &=\int d\mathbf {r} e^{-\gamma \mathbf {r} ^{\mathsf {T}}ww^{\mathsf {T}}\mathbf {r} }e^{-\mathbf {r} ^{\mathsf {T}}A'\mathbf {r} +\mathbf {s} '^{\mathsf {T}}\mathbf {r} }\left(w^{\mathsf {T}}\mathbf {r} \times w^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {r} ^{\mathsf {T}}}}\right)e^{-\mathbf {r} ^{\mathsf {T}}A\mathbf {r} +\mathbf {s} ^{\mathsf {T}}\mathbf {r} }\\&=\left(w^{\mathsf {T}}{\frac {\partial }{\partial \mathbf {s} '^{\mathsf {T}}}}\right)\times \left(w^{\mathsf {T}}\mathbf {s} -2w^{\mathsf {T}}A{\frac {\partial }{\partial \mathbf {s} ^{\mathsf {T}}}}\right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}\\&=\left({\frac {1}{2}}w^{\mathsf {T}}B'^{-1}\mathbf {v} \right)\times \left(w^{\mathsf {T}}\mathbf {s} -w^{\mathsf {T}}AB'^{-1}\mathbf {v} \right)e^{{\frac {1}{4}}\mathbf {v} ^{\mathsf {T}}B'^{-1}\mathbf {v} }\left({\frac {\pi ^{N}}{\det(B')}}\right)^{3/2}\,.\end{aligned}}$

where $B=A'+A$, $\mathbf {v} =\mathbf {s} '+\mathbf {s}$, $B'=B+\gamma ww^{\mathsf {T}}$.

Mathematical formulation for the ground state problem

Let us consider a time-independent physical system whose Hamiltonian H is Hermitian and bounded from below. We want to approximate the discrete eigenvalues of H and its wave functions

$H\Psi _{n}=E_{n}\Psi _{n}\;,$

where we ordered eigenvalues s.t.$E_{0}<E_{1}<...\;.$

It means that we would like to find such square integrable functions ${f_{i}}$, that $\langle (H-E_{i})f_{i}|(H-E_{i})f_{i}\rangle \ll {\mathit {eps}}\times E_{i}^{2}\langle f_{i}|f_{i}\rangle$, with some ${\mathit {eps}}\in \mathbf {R}$. Unfortunately in practice we don't know exact eigenvalues of the Hamiltonian, so first we have to find approximation to the energy $E_{i}$. The following theorem gives us the receipt. Here we would like to restrict ourselves to the ground state, but using the Min-max theorem one can extend to the whole discrete spectrum of the Hamiltonian

Theorem

The expectation value of the Hamiltonian for any $\phi$ from the state space is equal or larger then the ground state energy $E_{0}$.

Proof

Apparently the function $\phi$ can be decomposed in the orthogonal basis ${\Psi _{n}}$: $\phi =\sum a_{i}\Psi _{i}$. With this decomposition we write the mean value of the Hamiltonian: $E[\phi ]={\frac {<\phi |H|\phi >}{<\phi |\phi >}}={\frac {\sum _{i}|a_{i}|^{2}E_{i}}{\sum _{i}|a_{i}|^{2}}}$, from which follows that $E[\phi ]=E_{0}+{\frac {\sum _{i}|a_{i}|^{2}(E_{i}-E_{0})}{\sum _{i}|a_{i}|^{2}}}\geq E_{0}\;\Box$.

This statement is often called Ritz theorem and might be seen as a corollary of the Min-max theorem.

This result allows us compute an upper bound for the ground state energy.

The following theorem according to Weinstein allows us to rewrite our initial demand that $\langle (H-E_{i})f_{i}|(H-E_{i})f_{i}\rangle \ll {\mathit {eps}}\times E_{i}^{2}\langle f_{i}|f_{i}\rangle$ in terms of the variance $\sigma ^{2}[\phi ]={\frac {<\phi |(H-E)^{2}|\phi >}{<\phi |\phi >}}$

Theorem

There exist at least one eigenvalue in the interval ${\bigg [}E[\phi ]-\sigma [\phi ],E[\phi ]+\sigma [\phi ]{\bigg ]}$ .

Proof

We write $\phi$ in the ${\Psi _{n}}$ basis, and get $\sigma ^{2}[\phi ]={\frac {\sum _{i}|a_{i}|^{2}(E_{i}-E)^{2}}{\sum _{i}|a_{i}|^{2}}}$. There exist integer $k$, s.t. $(E_{k}-E)^{2}<(E_{i}-E)^{2},\forall i\neq k$. With this we rewrite variance $\sigma ^{2}[\phi ]=(E_{k}-E)^{2}+{\frac {\sum _{i}|a_{i}|^{2}{\bigg (}(E_{i}-E)^{2}-(E_{k}-E)^{2}{\bigg )}}{\sum _{i}|a_{i}|^{2}}}\geq (E_{k}-E)^{2}\;\Box .$

This result might be useful if and only if the lower bound can be calculated as close as possible to the ground state energy.

With these theorems we see the way to proceed:

1. Take convenient basis in the state space of the Hamiltonian.

2. Cut the basis size to some finite number.

3. Minimise expectation value of the Hamiltonian in this basis.

4. Enlarge basis and do step 3.

5. Do steps 3,4 as long as needed to insure convergence of the ground state energy.

6. Calculate variance.

7. If variance is larger than some precision value than enlarge basis size and do 3,4,5,6 again, otherwise we are done.

In practice steps 3,4,5 alone can give accurate value of energy. Steps 6,7 are needed for approximation of the wave function. This is due to the following theorem

Theorem

The expectation value of the Hamiltonian is stationary in the neighbourhood of the discrete eigenvalues.

Proof

So in general it is easier to get accurate approximation to the energy than to other observables.

Basis

We want to start with the first step: take some convenient basis. We would like to define convenient for our problem

1. Simple transformation from one system of coordinates to another.

2. Possibility to eliminate the centre of mass.

3. Easy computations for the overlap and kinetic energy.

It is of advantage to introduce rescaled coordinates,

$\mathbf {q} _{i}={\sqrt {\frac {m_{i}}{m}}}\mathbf {r} _{i}\;,$

where $m$ is a conveniently chosen mass scale. Indeed the kinetic energy $T$ and the harmonic trap $V_{h}$ have a more symmetric form in the rescaled coordinates,

$T=-{\frac {\hbar ^{2}}{2m}}\sum _{i}{\frac {\partial }{\partial \mathbf {q} _{i}^{2}}}\;,V_{h}={\frac {m\omega ^{2}}{2}}\sum _{i=1}^{N}\mathbf {q} _{i}^{2}\;.$

The Jacobian of the transformation from $\mathbf {r}$ to $\mathbf {q}$ is equal

${\frac {\partial (\mathbf {q} _{1},\dots ,\mathbf {q} _{N})}{\partial (\mathbf {r} _{1},\dots ,\mathbf {r} _{N})}}=\prod _{i=1}^{N}\left({\frac {m_{i}}{m}}\right)^{3/2}\;.$

A further suitable linear transformation to a new set of coordinates is possible,

$\mathbf {x} _{i}=\sum _{j=1}^{N}U_{ij}\mathbf {q} _{j}\;,$

or, in matrix notation,

$\mathbf {x} =U\mathbf {q} \;,$

where $U$ is the transformation matrix.

If the transformation matrix is unitary, $U^{T}U=1$, the diagonal form of the kinetic energy and the harmonic trap is preserved in the new coordinates,

$T=-{\frac {\hbar ^{2}}{2m}}\sum _{i}{\frac {\partial }{\partial \mathbf {x} _{i}^{2}}}\;,V_{h}={\frac {m\omega ^{2}}{2}}\sum _{i=1}^{N}\mathbf {x} _{i}^{2}\;.$

Last transformation is of particular use if new system has the coordinate

$\mathbf {x_{N}} \sim \sum _{i}\mathbf {r_{i}} \;,$

which can be seen as a center of mass coordinate. It allows us to work with a wave function in the form

$\Psi (\mathbf {r_{1}} ,\mathbf {r_{2}} ,...,\mathbf {r_{N}} )=\Phi (\mathbf {x_{1}} ,\mathbf {x_{2}} ,...,\mathbf {x_{N-1}} )\phi (\mathbf {x_{N}} )\;,$

where $\phi (\mathbf {x_{N}} )$ is the ground state wave function for the oscillator potential.

Correlated Gaussians General Case

First we consider trial wave function in the basis of completely general shifted Gaussian, which can be used to describe a system in the external field with anisotropic inter-particle interaction

$\left|\phi \right\rangle =\sum _{k=1}^{K}C_{k}\left|A^{(k)},s^{(k)};x\right\rangle$

where

$\left|A,s;x\right\rangle =\exp \left(-\sum _{i,j=1}^{D*n}A_{ij}(x_{i}-s_{i})(x_{j}-s_{j})\right)\equiv e^{-(x-s)^{T}A(x-s)}\;,$

$A$, a symmetric positive-defined matrix, and $s$, a shift vector, are the non-linear parameters of the Gaussian and n=N-1. With this definition we have $K({\frac {1}{2}}D*n(D*n+1)+D*n)$ non-linear variational parameters. To find those one can use deterministic methods (e.g. Powell's method) or methods based on a stochastic search. We use latter approach though we find $K$ linear variational parameters through a full minimization with respect to a given set of non-linear parameters.

Matrix elements

The matix elements can be determined analytically either by diagonalizing the matrix $B$ or Cholesky decompose $B=LL^{T}$ (since it is positive definite) and change basis to coordinates where the matrix is either diagonal or unit. This way many integrals can be determined by iterative integration.

Overlap

Correlated Gaussians are generally non-orthogonal and the overlap is therefore non-diagonal,

$N\equiv \langle A',s';x|A,s;x\rangle =\int d^{D}x_{1}\dots \int d^{D}x_{n}e^{-(x-s)^{T}A(x-s)-(x-s')^{T}A'(x-s')}={\frac {\pi ^{\frac {D*n}{2}}}{\sqrt {\det B}}}e^{(-s^{T}As-s'^{T}A's'+{\frac {1}{4}}v^{T}B^{-1}v)}\;.$

where we defined $B=A+A'\;,v=2As+2A's'$

Kinetic energy

Here we calculate kinetic energy

$-\langle A',s';x|{\frac {\partial }{\partial x}}^{T}\Lambda {\frac {\partial }{\partial x}}|A,s;x\rangle =N\times \left(2{\rm {tr}}(A'\Lambda AB^{-1})+4u^{T}A'\Lambda Au-4u^{T}(A'\Lambda As+A\Lambda A's')+4s'^{T}A'\Lambda As\right)$

where we defined $u={\frac {1}{2}}B^{-1}v$ for $\Lambda ={\frac {1}{2}}I$ with $I$ to be the identity matrix, one can get simpler expression after noticing that

$-{\frac {1}{2}}{\frac {\partial }{\partial x}}^{T}{\frac {\partial }{\partial x}}|A,s;x\rangle =({\rm {tr}}A-2(x-s)^{T}AA(x-s))|A,s;x\rangle$

To proceed one has to derive the following identities:

$\langle A',s';x|\left(a^{T}x\right)|A,s;x\rangle =N\times (a^{T}u)\;,$

$\langle A',s';x|\left(x^{T}Fx\right)|A,s;x\rangle =N\times \left(u^{T}Fu+{\frac {1}{2}}{\rm {tr}}(FB^{-1})\right)\;.$

To calculate

$\langle A',s';x|(-{\frac {1}{2}}{\frac {\partial }{\partial x}}^{T}{\frac {\partial }{\partial x}})(-{\frac {1}{2}}{\frac {\partial }{\partial x}}^{T}{\frac {\partial }{\partial x}})|A,s;x\rangle$

which we need to calculate variance we have to calculate the following matrix element

${\begin{array}{lcl}\langle A',s';x|(x^{T}Dx)(x^{T}D_{1}x)|A,s;x\rangle =N\times &{\bigg (}[u^{T}Du][u^{T}D_{1}u]+{\frac {1}{2}}[u^{T}D_{1}u]{\rm {tr}}(DB^{-1})+{\frac {1}{2}}[u^{T}Du]{\rm {tr}}(D_{1}B^{-1})+\\&2u^{T}DB^{-1}D_{1}u+{\frac {1}{4}}{\rm {tr}}(DB^{-1}){\rm {tr}}(D_{1}B^{-1})+{\rm {tr}}(DB^{-1}D_{1}B^{-1}){\bigg )}\;,\end{array}}$

Potential energy

Here we calculate matrix element

$\langle A',s';x|V|A,s;x\rangle =\sum _{i<j}^{N}\int d^{D}x_{1}\dots \int d^{D}x_{n}V_{ij}e^{-(x-s)^{T}A(x-s)-(x-s')^{T}A'(x-s')}\;.$

In general we can not write analytical expression for this integral, but we can reduce it to D dimensional integrals. For example, consider just one term from the sum

$I=\int d^{D}x_{1}\dots \int d^{D}x_{n}V_{ij}e^{-(x-s)^{T}A(x-s)-(x-s')^{T}A'(x-s')}\;,$

to simplify this integral we have to make transformation from Jacobi set $(\mathbf {x_{1}} ,\mathbf {x_{2}} ,...,\mathbf {x_{n}} )$ , where matrices A,A',s,s' are defined to the Jacobi set $(\mathbf {y_{1}} ,\mathbf {y_{2}} ,...,\mathbf {y_{n}} )$, where $\mathbf {y_{1}} =\mathbf {r_{i}} -\mathbf {r_{j}}$. transformation between those sets are provided through the orthogonal matrix U: x=Uy. With this we write

$I=\int d^{D}y_{1}\dots \int d^{D}y_{n}V_{ij}(\mathbf {y_{1}} )e^{-(y-U^{T}s)^{T}U^{T}AU(y-U^{T}s)-(y-U^{T}s')^{T}U^{T}A'U(y-U^{T}s')}=\int d^{D}y_{1}V_{ij}(\mathbf {y_{1}} )f(\mathbf {y_{1}} )\;,$

where

$f(\mathbf {y_{1}} )=\int d^{D}y_{2}\dots \int d^{D}y_{n}e^{-(y-U^{T}s)^{T}U^{T}AU(y-U^{T}s)-(y-U^{T}s')^{T}U^{T}A'U(y-U^{T}s')}\;,$

can be found analytically.

If we can write potential as a sum of Gaussians $V_{ij}(\mathbf {y_{1}} )=\sum _{k}c_{k}e^{-\sum _{a,b=1}^{a,b=D}(y_{a}-u_{a}^{k})^{T}F_{ab}^{k}(y_{b}-u_{b}^{k})}$ then the integral $I$ can be found analytically in the same way as we found overlap.

Particles with spin

To consider particles with spin we add spin part to the trial wave function

$\left|\phi \right\rangle =\sum _{k=1}^{K}C_{k}\left|A^{(k)},s^{(k)};x\right\rangle \chi _{k}\;$

where for particles with spin = 1/2 function $\chi _{k}$ is just an array of N elements. Each element is an eigenfunction of the spin's projection on the predefined axis. For example $\chi _{k}=|{\frac {1}{2}},{\frac {1}{2}},...,{\frac {1}{2}}>$ defines the system with all particles have spin in the same direction. Next we define the spin operator $\mathbf {s^{(i)}}$ that acts on the particle number $i$ in the following way

$<\pm {\frac {1}{2}}|s_{z}^{(i)}|\pm {\frac {1}{2}}>=\pm {\frac {1}{2}}$

$<\pm {\frac {1}{2}}|s_{x}^{(i)}|\mp {\frac {1}{2}}>={\frac {1}{2}}$

$<\pm {\frac {1}{2}}|s_{y}^{(i)}|\mp {\frac {1}{2}}>=\mp {\frac {\mathbf {i} }{2}}$

and zero otherwise.

Spin-orbit

Here we discuss spin-orbit potential of the form $V_{ij,LS}=V(|\mathbf {r_{i}} -\mathbf {r_{j}} |)\mathbf {L^{(ij)}} {\mathbf {S^{(ij)}} }$, where $\mathbf {S^{(ij)}} =\mathbf {s^{(i)}} +\mathbf {s^{(j)}}$ and $L_{k}^{(ij)}=-\mathbf {i} \sum _{l,m=1}^{l,m=D}e_{klm}y_{l}{\frac {\partial }{\partial y_{m}}}$ - relative angular momentum, where $e_{klm}$ - Levi-Chevita symbol, and $\mathbf {y} =\mathbf {r_{i}} -\mathbf {r_{j}}$. We have to calculate following matrix element

$I_{k}=\left\langle A',s';x|V(y)e_{klm}y_{l}{\frac {\partial }{\partial y_{m}}}|A,s;x\right\rangle$

again we are making transformation from the Jacobi set $(\mathbf {x_{1}} ,...,\mathbf {x_{n}} )$ to the Jacobi set $(\mathbf {y_{1}} =\mathbf {y} ,...,\mathbf {y_{n}} )$ using a transformation matrix U: x=Uy.

${\begin{array}{lcl}I_{k}=&\int d^{D}y_{1}\dots \int d^{D}y_{n}V(\mathbf {y_{1}} )e_{klm}y_{l}{\bigg (}-2(U^{T}AU)_{jm}(y-U^{T}s)_{j}{\bigg )}\\&e^{-(y-U^{T}s)^{T}U^{T}AU(y-U^{T}s)-(y-U^{T}s')^{T}U^{T}A'U(y-U^{T}s')},0<j\leq n*D,0<k,l,m\leq D\;,\end{array}}$

Correlated Gaussians with Super Vectors

In the previous section we considered completely general setup, which is suitable for any inter-particle potentials and external fields. This approach is far from optimal if for example we are interested in ground state of N bosons with isotropic pairwise interaction, because in this case we know that our ground state must have zero orbital momentum, with this in mind we write the trial wave function in the smaller variational basis:

$\left|\phi \right\rangle =\sum _{k=1}^{K}C_{k}\left|A^{(k)},\mathbf {s^{(k)}} ;x\right\rangle$

where

$\left|A,\mathbf {s} ;x\right\rangle =\exp \left(-\sum _{i,j=1}^{n}A_{ij}(\mathbf {x_{i}} -\mathbf {s_{i}} )\cdot (\mathbf {x_{j}} -\mathbf {s_{j}} )\right)\equiv e^{-(\mathbf {x} -\mathbf {s} )^{T}A(\mathbf {x} -\mathbf {s} )}\;.$

If we put shift vectors to be zeros $\mathbf {s_{i}} =0,\mathbf {s_{j}} =0$, then the trial wave function treats Cartesian components of vectors $\mathbf {x}$ equivalently, which leads to zero angular momentum, otherwise the wave function will contain all possible angular momentum and we need an effective procedure to build an eigenstate for a given angular momentum. Matrix elements for this trial wave functions can be obtained from the general case, but we write it explicitly.

Matrix Elements

Overlap

$N\equiv \langle A',\mathbf {s'} ;x|A,\mathbf {s} ;x\rangle =\int d^{D}x_{1}\dots \int d^{D}x_{n}e^{-(\mathbf {x} -\mathbf {s} )^{T}A(\mathbf {x} -\mathbf {s} )-(\mathbf {x} -\mathbf {s'} )^{T}A'(\mathbf {x} -\mathbf {s'} )}={\bigg (}{\frac {\pi ^{n}}{\det B}}{\bigg )}^{\frac {D}{2}}e^{(-\mathbf {s} ^{T}A\mathbf {s} -\mathbf {s'} ^{T}A'\mathbf {s'} +{\frac {1}{4}}\mathbf {v} ^{T}B^{-1}\mathbf {v} )}\;.$

where we defined $B=A+A'\;,\mathbf {v} =2A\mathbf {s} +2A'\mathbf {s'}$

Kinetic energy

${\begin{array}{lcl}-\langle A',\mathbf {s'} ;x|{\frac {\partial }{\partial x}}^{T}\Lambda {\frac {\partial }{\partial x}}|A,\mathbf {s} ;x\rangle =&N\times \left(2D{\rm {tr}}(A'\Lambda AB^{-1})+4\mathbf {u} ^{T}A'\Lambda A\mathbf {u} -4\mathbf {u} ^{T}(A'\Lambda A\mathbf {s} +A\Lambda A'\mathbf {s'} )+4\mathbf {s'} ^{T}A'\Lambda A\mathbf {s} \right)=\\&N\times \left(2D{\rm {tr}}(A'\Lambda AB^{-1})+4(\mathbf {u} -\mathbf {s'} )^{T}A'\Lambda A(\mathbf {u} -\mathbf {s} )\right)\end{array}}$

where we defined $\mathbf {u} ={\frac {1}{2}}B^{-1}\mathbf {v}$

Angular momentum

We consider the matrix element of the operator $O=\sum _{i,j=1}^{i,j=n}a_{i}b_{j}[\mathbf {x_{i}} \times \mathbf {\frac {\partial }{\partial \mathbf {x_{j}} }} ]\equiv [a\mathbf {x} \times b{\frac {\partial }{\partial \mathbf {x} }}]$, choice of $a,b$ can give the total angular momentum or an angular momentum for appropriate relative coordinate.

First we calculate matrix elements of the form

$\langle A',\mathbf {s'} ;x|a\mathbf {x} |A,\mathbf {s} ;x\rangle =N\cdot (a\mathbf {u} )$

$\langle A',\mathbf {s'} ;x|[a\mathbf {x} \times b\mathbf {x} ]|A,\mathbf {s} ;x\rangle =N\cdot (a\mathbf {u} \times b\mathbf {u} )$

and now we can calculate matrix element for operator $O$

$\langle A',\mathbf {s'} ;x|[a\mathbf {x} \times b{\frac {\partial }{\partial \mathbf {x} }}]|A,\mathbf {s} ;x\rangle =-2\langle A',\mathbf {s'} ;x|[a\mathbf {x} \times bA\mathbf {x} ]|A,\mathbf {s} ;x\rangle +2\langle A',\mathbf {s'} ;x|[a\mathbf {x} \times bA\mathbf {s} ]|A,\mathbf {s} ;x\rangle =2N[a\mathbf {u} \times bA(\mathbf {s} -\mathbf {u} )]\;.$

We define total angular momentum to be $\mathbf {L} \equiv \sum _{i}^{N}(\mathbf {r} _{i}\times \mathbf {p} _{i})$. If we make transformation to the Jacobi set, than we obtain $\mathbf {L} =\sum _{i}^{N}(\mathbf {x} _{i}\times \mathbf {P} _{i})$, where $\mathbf {P} _{i}$ is the linear momentum corresponding to the $i$ coordinate. So if we assume that the system as a whole is at rest s.t. $\mathbf {P} _{N}=0$ than the following matrix element defines total angular momentum

$<L^{2}>=-\hbar ^{2}\sum _{i=1}^{n}\sum _{k=1}^{n}\langle A',\mathbf {s'} ;x|(\mathbf {x} _{i}\times {\frac {\partial }{\partial \mathbf {x} _{i}}})(\mathbf {x} _{k}\times {\frac {\partial }{\partial \mathbf {x} _{k}}})|A,\mathbf {s} ;x\rangle =\hbar ^{2}\sum _{i=1}^{n}\sum _{k=1}^{n}\langle A',\mathbf {s'} ;x|\leftarrow (\mathbf {x} _{i}\times {\frac {\partial }{\partial \mathbf {x} _{i}}})(\mathbf {x} _{k}\times {\frac {\partial }{\partial \mathbf {x} _{k}}})|A,\mathbf {s} ;x\rangle$

After simplification (first we rotate to the set of coordinate, where matrix $A$ takes a diagonal form, than we rotate back and rotate to the set where matrix $A'$ takes a diagonal form, and again rotate back) we obtain

$<L^{2}>=4\hbar ^{2}\sum _{i,m,k,l=1}^{n}\langle A',\mathbf {s'} ;x|A_{im}'(\mathbf {x} _{i}\times \mathbf {s'} _{m})A_{kl}(\mathbf {x} _{k}\times \mathbf {s} _{l})|A,\mathbf {s} ;x\rangle$

We take the following integral

$\langle A',\mathbf {s'} ;x|x_{m}^{l}x_{k}^{f}|A,\mathbf {s} ;x\rangle =N{\bigg (}{\frac {1}{2}}(B^{-1})_{km}\delta ^{fl}+{\frac {1}{4}}(B^{-1})_{kc}v_{c}^{f}(B^{-1})_{mn}v_{n}^{l}{\bigg )};l,f=1,2,3;m,k=1,2,...,n$

where $\delta ^{fl}$ - is a Kronecker delta.

With this we write total angular momentum

$<L^{2}>=4\hbar ^{2}N{\bigg (}\mathbf {s'} ^{T}A'B^{-1}A\mathbf {s} +{\frac {1}{4}}\sum _{k,c}(A'B^{-1})_{kc}[\mathbf {v} _{k}\times \mathbf {s'} _{c}]\cdot \sum _{a,b}(AB^{-1})_{ab}[\mathbf {v} _{a}\times \mathbf {s} _{b}]{\bigg )}$

Appendix A: some integrals with 3D Gaussians

Not-shifted

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}=\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}{\vec {x}}^{2}={\frac {3}{2}}{\frac {1}{\beta }}\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}({\vec {a}}{\vec {x}})({\vec {b}}{\vec {x}})=({\vec {a}}{\vec {b}}){\frac {1}{3}}\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}{\vec {x}}^{2}=({\vec {a}}{\vec {b}}){\frac {1}{2\beta }}\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}{\frac {({\vec {a}}{\vec {x}})({\vec {b}}{\vec {x}})}{{\vec {x}}^{2}}}=({\vec {a}}{\vec {b}}){\frac {1}{3}}\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}=({\vec {a}}{\vec {b}}){\frac {1}{3}}\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta x^{2}}\left(3{\frac {({\vec {a}}{\vec {x}})({\vec {b}}{\vec {x}})}{{\vec {x}}^{2}}}-({\vec {a}}{\vec {b}})\right)=0$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}({\vec {x}}^{2})^{2}={\frac {15}{4}}{\frac {\pi ^{3/2}}{\beta ^{7/2}}}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}}({\vec {x}}^{2})^{3}={\frac {105}{8}}{\frac {\pi ^{3/2}}{\beta ^{9/2}}}$

---

$\int _{-\infty }^{+\infty }d^{3}xe^{-\alpha ^{\prime }{\vec {x}}^{2}}e^{-\alpha {\vec {x}}^{2}}=\left({\frac {\pi }{\alpha +\alpha ^{\prime }}}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\alpha ^{\prime }{\vec {x}}^{2}}{\vec {x}}^{2}e^{-\alpha {\vec {x}}^{2}}={\frac {3}{2}}{\frac {1}{\alpha +\alpha ^{\prime }}}\left({\frac {\pi }{\alpha +\alpha ^{\prime }}}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\alpha ^{\prime }{\vec {x}}^{2}}\left(-{\frac {\partial ^{2}}{\partial {\vec {x}}^{2}}}\right)e^{-\alpha {\vec {x}}^{2}}=6{\frac {\alpha \alpha ^{\prime }}{\alpha +\alpha ^{\prime }}}\left({\frac {\pi }{\alpha +\alpha ^{\prime }}}\right)^{3/2}$

---

n-body

$\int _{-\infty }^{+\infty }d^{3}x_{1}\dots d^{3}x_{n}\exp \left(-\sum _{i=1}^{n}\alpha _{i}^{\prime }{\vec {x}}_{i}^{2}\right)\exp \left(-\sum _{i=1}^{n}\alpha _{i}{\vec {x}}_{i}^{2}\right)=\left({\frac {\pi ^{n}}{\prod _{i=1}^{n}(\alpha _{i}+\alpha _{i}^{\prime })}}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}x_{1}\dots d^{3}x_{n}\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}^{\prime }{\vec {x}}_{j}\right)\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}{\vec {x}}_{j}\right)=\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2},B=A+A^{\prime }$

$\int _{-\infty }^{+\infty }d^{3}x_{1}\dots d^{3}x_{n}\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}^{\prime }{\vec {x}}_{j}\right)\sum _{ij}{\vec {x}}_{i}F_{ij}{\vec {x}}_{j}\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}{\vec {x}}_{j}\right)=\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}{\frac {3}{2}}\;\mathrm {trace} (FB^{-1}),B=A+A^{\prime }$

$\int _{-\infty }^{+\infty }d^{3}x_{1}\dots d^{3}x_{n}\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}^{\prime }{\vec {x}}_{j}\right)\sum _{ij}\left(-{\frac {\partial }{\partial {\vec {x}}_{i}}}\Lambda _{ij}{\frac {\partial }{\partial {\vec {x}}_{j}}}\right)\exp \left(-\sum _{ij}{\vec {x}}_{i}A_{ij}{\vec {x}}_{j}\right)=\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}6\;\mathrm {trace} (A^{\prime }\Lambda AB^{-1}),B=A+A^{\prime }$

Shifted

${\vec {u}}\doteq {\frac {1}{2}}\beta ^{-1}{\vec {v}}\;,$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}+{\vec {v}}{\vec {x}}}=e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}+{\vec {v}}{\vec {x}}}{\vec {a}}{\vec {x}}=\left({\vec {a}}{\frac {\partial }{\partial {\vec {v}}}}\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}={\frac {1}{2\beta }}\left({\vec {a}}{\vec {v}}\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}=\left({\vec {a}}{\vec {u}}\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}\;,$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}+{\vec {v}}{\vec {x}}}({\vec {a}}{\vec {x}})({\vec {b}}{\vec {x}})=\left({\vec {a}}{\frac {\partial }{\partial {\vec {v}}}}\right)\left({\vec {b}}{\frac {\partial }{\partial {\vec {v}}}}\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}=\left({\frac {1}{2\beta }}({\vec {a}}{\vec {b}})+{\frac {1}{4\beta ^{2}}}({\vec {a}}{\vec {v}})({\vec {b}}{\vec {v}})\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d^{3}xe^{-\beta {\vec {x}}^{2}+{\vec {v}}{\vec {x}}}{\vec {x}}^{2}=\left({\frac {3}{2}}{\frac {1}{\beta }}+{\vec {u}}^{2}\right)e^{{\frac {1}{4}}{\frac {{\vec {v}}^{2}}{\beta }}}\left({\frac {\pi }{\beta }}\right)^{3/2}$

---

Here we need to indroduce a notation to simplify the expressions:

${\tilde {x}}\doteq \left({\begin{array}{c}{\vec {x}}_{1}\\\vdots \\{\vec {x}}_{n}\end{array}}\right)\;,\;({\tilde {x}})_{i}={\vec {x}}_{i}\;,\;\left(B{\tilde {x}}\right)_{i}\doteq \sum _{j}B_{ij}{\vec {x}}_{j}\;,\;\left({\tilde {x}}^{T}B\right)_{i}\doteq \sum _{j}{\vec {x}}_{j}B_{ji}\;,\;{\tilde {v}}^{T}{\tilde {x}}\doteq \sum _{i=1}^{n}{\vec {v}}_{i}\cdot {\vec {x}}_{i}\;.$

$d{\tilde {x}}\doteq d^{3}x_{1}\dots d^{3}x_{n}\;.$

---

$\int _{-\infty }^{+\infty }d{\tilde {x}}e^{-{\tilde {x}}^{T}B{\tilde {x}}+{\tilde {v}}^{T}{\tilde {x}}}=e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}$

$\int _{-\infty }^{+\infty }d{\tilde {x}}e^{-{\tilde {x}}^{T}B{\tilde {x}}+{\tilde {v}}^{T}{\tilde {x}}}\left({\tilde {a}}^{T}{\tilde {x}}\right)=\left({\tilde {a}}^{T}{\frac {\partial }{\partial {\tilde {v}}^{T}}}\right)e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}=\left({\tilde {a}}^{T}{\tilde {u}}\right)e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}\;,\;{\tilde {u}}={\frac {1}{2}}B^{-1}{\tilde {v}}$

$\int _{-\infty }^{+\infty }d{\tilde {x}}e^{-{\tilde {x}}^{T}B{\tilde {x}}+{\tilde {v}}^{T}{\tilde {x}}}\left({\tilde {a}}^{T}{\tilde {x}}\right)\left({\tilde {b}}^{T}{\tilde {x}}\right)=\left({\tilde {a}}^{T}{\frac {\partial }{\partial {\tilde {v}}^{T}}}\right)\left({\tilde {b}}^{T}{\frac {\partial }{\partial {\tilde {v}}^{T}}}\right)e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}=\left({\frac {1}{2}}{\tilde {a}}^{T}B^{-1}{\tilde {b}}+({\tilde {a}}^{T}{\tilde {u}})({\tilde {b}}^{T}{\tilde {u}})\right)e^{{\frac {1}{2}}{\tilde {v}}^{T}{\tilde {u}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}\;,\mathrm {where} \;B=B^{T}$

$\int _{-\infty }^{+\infty }d{\tilde {x}}e^{-{\tilde {x}}^{T}B{\tilde {x}}+{\tilde {v}}^{T}{\tilde {x}}}\left({\tilde {x}}^{T}F{\tilde {x}}\right)=\left({\frac {3}{2}}\mathrm {trace} (FB^{-1})+{\tilde {u}}^{T}F{\tilde {u}}\right)e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}\;,\;$

${\begin{aligned}&\int _{-\infty }^{+\infty }d{\tilde {x}}\;e^{-{\tilde {x}}^{T}A'{\tilde {x}}+{\tilde {s}}'^{T}{\tilde {x}}}\left(-{\frac {\partial }{\partial {\tilde {x}}}}\Lambda {\frac {\partial }{\partial {\tilde {x}}^{T}}}\right)e^{-{\tilde {x}}^{T}A{\tilde {x}}+{\tilde {s}}^{T}{\tilde {x}}}\\&=\left(6\mathrm {trace} (A'\Lambda AB^{-1})+4{\tilde {u}}^{T}A'\Lambda A{\tilde {u}}-2{\tilde {u}}^{T}A'\Lambda {\tilde {s}}-2{\tilde {s}}'^{T}\Lambda A{\tilde {u}}+{\tilde {s}}'^{T}\Lambda {\tilde {s}}\right)e^{{\frac {1}{4}}{\tilde {v}}^{T}B^{-1}{\tilde {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)^{3/2}\;,\;B=A+A'\;,\;{\tilde {v}}={\tilde {s}}+{\tilde {s}}'\end{aligned}}$

Appendix B: some integrals with 1D Gaussians

First the overlap (denoted N and used in the following expressions)

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}=e^{{\frac {1}{4}}{\vec {v}}B^{-1}{\vec {v}}}\left({\frac {\pi ^{n}}{\det(B)}}\right)=N\end{aligned}}$

Polynomial like terms

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times \left({\vec {a}}^{T}{\vec {x}}\right)=N\times {\frac {1}{2}}\left({\vec {a}}^{T}B^{-1}{\vec {v}}\right)\end{aligned}}$

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times \left({\vec {x}}^{T}F{\vec {x}}\right)=N\times \left[{\frac {1}{2}}Trace\left(FB^{-1}\right)+{\frac {1}{4}}{\vec {v}}^{T}B^{-1}{\vec {v}}\right]\end{aligned}}$

for any matrix F with the right dimensions including $F={\vec {a}}{\vec {b}}^{T}$.

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times \left({\vec {a}}^{T}{\vec {x}}{\vec {x}}^{T}F{\vec {x}}\right)=\end{aligned}}$

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times \left({\vec {x}}^{T}F{\vec {x}}\right)\left({\vec {x}}^{T}P{\vec {x}}\right)=\end{aligned}}$

Exponential and sinusoidal terms

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times ae^{{\vec {c}}^{T}{\vec {x}}}=N\times ae^{{\frac {1}{4}}{\vec {c}}^{T}B^{-1}{\vec {c}}+{\frac {1}{2}}{\vec {c}}B^{-1}{\vec {v}}}\end{aligned}}$

(Co)Sine term:

Due to the Euler formulae for (Co)Sine and the linearity we can express these matrix elements as the (Real) Imaginary part of the above matrix element with complex vector ${\vec {c}}=i{\vec {l}}$

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times a\cos \left({\vec {l}}^{T}{\vec {x}}\right)=N\times a\cos \left({\frac {1}{2}}{\vec {l}}^{T}B^{-1}{\vec {v}}\right)e^{{\frac {-1}{4}}{\vec {l}}B^{-1}{\vec {l}}}\end{aligned}}$

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times a\sin \left({\vec {l}}^{T}{\vec {x}}\right)=N\times a\sin \left({\frac {1}{2}}{\vec {l}}^{T}B^{-1}{\vec {v}}\right)e^{{\frac {-1}{4}}{\vec {l}}B^{-1}{\vec {l}}}\end{aligned}}$

---

Gaussian term:

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times G_{0}{\sqrt {\frac {g}{\pi }}}e^{-g\left({\vec {a}}^{T}{\vec {x}}-R\right)^{2}}=N\times {\sqrt {\frac {g}{\pi \left(1+g{\vec {a}}^{T}B^{-1}{\vec {a}}\right)}}}e^{{\frac {-g}{1+g{\vec {a}}^{T}B^{-1}{\vec {a}}}}\left({\frac {1}{2}}{\vec {a}}^{T}B^{-1}{\vec {v}}-R\right)^{2}}\end{aligned}}$

Delta term:

This can be calculated as the limiting case $g\rightarrow \infty$ of the previous matrix element

${\begin{aligned}\int _{-\infty }^{\infty }\,d{\vec {x}}e^{-{\vec {x}}^{T}B{\vec {x}}+{\vec {v}}^{T}{\vec {x}}}\times \delta \left({\vec {a}}^{T}{\vec {x}}-R\right)=N\times {\sqrt {\frac {1}{\pi {\vec {a}}^{T}B^{-1}{\vec {a}}}}}e^{{\frac {-1}{{\vec {a}}^{T}B^{-1}{\vec {a}}}}\left({\frac {1}{2}}{\vec {a}}^{T}B^{-1}{\vec {v}}-R\right)^{2}}\end{aligned}}$