# Analytic matrix elements with shifted correlated Gaussians

D.V. Fedorov
Institute of Physics and Astronomy, Aarhus University,
Ny Munkegade 120, 8000 Aarhus C, Denmark
<EMAIL_ADDRESS>
###### Abstract

Matrix elements between shifted correlated Gaussians of various potentials
with several form-factors are calculated analytically. Analytic matrix
elements are of importance for the correlated Gaussian method in quantum few-
body physics.

##  1 Introduction

Correlated Gaussian method is a popular variational method to solve quantum-
mechanical few-body problems in molecular, atomic, and nuclear physics [1, 2].
One of the important advantages of the correlated Gaussian method is the ease
of computing the matrix elements. In some cases the matrix elements and even
their gradients with respect to optimization parameters are fully analytic [2,
3, 4, 5]. This enables extensive numerical optimizations to be carried out
leading to accurate results [6, 7, 8] despite the incorrect functional form of
the Gaussians in certain asymptotic regions of the configuration space.

Although a number of analytic matrix elements have been calculated for
different potentials and different forms of correlated Gaussians, one
combination — short-range potentials with shifted correlated Gaussians — is
still missing [1]. Indeed historically the shifted Gaussians have been applied
more often to Coulombic systems rather than to atomic and nuclear systems
where the short-range interactions are most important.

In this paper several types of short-range potentials are considered in a
search for the form-factors that produce analytic matrix elements with shifted
correlated Gaussians. A brief introduction to the correlated Gaussian method
is given first and then the analytic matrix elements are introduced.

##  2 Correlated Gaussian method

Correlated Gaussian method is a variational method where the coordinate part
of the wave-function of a quantum few-body system is expanded in terms of
correlated Gaussians. Various forms of correlated Gaussians have been
concieved [1], one of them being the shifted correlated Gaussian, $|g\rangle$,
which for a system of $N$ particles with coordinates $\vec{r}_{i}|_{i=1\dots
N}$ has the form

$\langle\mathbf{r}|g\rangle=\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)\equiv
e^{-\mathbf{r}^{\mathsf{T}}A\mathbf{r}+\mathbf{s}^{\mathsf{T}}\mathbf{r}}\,,$
(1)

where $\mathbf{r}$ is size-$N$ column of particle coordinates $\vec{r}_{i}$,
$A=\\{A_{ij}\\}$ is a size-$N$ square symmetric positive-definite correlation
matrix, $\mathbf{s}$ is size-$N$ column of shift vectors $\vec{s}_{i}$, and
where the following notation has been introduced,

$\mathbf{r}^{\mathsf{T}}A\mathbf{r}\equiv\sum_{i,j=1}^{N}A_{ij}\,\vec{r}_{i}\cdot\vec{r}_{j}\,,\quad\mathbf{s}^{\mathsf{T}}\mathbf{r}\equiv\sum_{i=1}^{N}\vec{s}_{i}\cdot\vec{r}_{i}\,,$
(2)

where “$\cdot$” denotes the dot-product of two vectors. The elements of the
correlation matrix and the shift vectors are the non-linear parameters of the
Gaussians.

The coordinate part of the few-body wave-function $|\psi\rangle$ is
represented as a linear combination of several Gaussians,

$|\psi\rangle=\sum_{i=1}^{n_{g}}c_{i}|g_{i}\rangle\,,$ (3)

where the coefficients $c_{i}$ are the linear parameters, and where $n_{g}$ is
the number of Gaussians. Inserting this representation into the Schrödinger
equation,

$\hat{H}|\psi\rangle=E|\psi\rangle\,,$ (4)

where $\hat{H}$ is the Hamiltonian of the few-body system, and multiplying
from the left with $\langle g_{j}|$ leads to the generalized matrix eigenvalue
equation,

$\mathcal{H}c=E\mathcal{N}c\,,$ (5)

where $c=\\{c_{i}\\}$ is the column of the linear parameters, and where the
Hamilton matrix $\mathcal{H}=\\{\mathcal{H}_{ij}\\}$ and the overlap matrix
$\mathcal{N}=\\{\mathcal{N}_{ij}\\}$ are given as

$\mathcal{H}_{ij}=\langle
g_{i}|\hat{H}|g_{j}\rangle\,,\quad\mathcal{N}_{ij}=\langle
g_{i}|g_{j}\rangle\,.$ (6)

The linear parameters together with the energy spectrum are found by solving
the generalized eigenvalue problem (5) numerically using the standard linear
algebra methods [9]. The non-linear parameters are optimized by one of many
optimization methods which typically involve elements of stochastic-
evolutionary [10] and direct optimization algorithms [7, 8]. The direct
optimization algorithms often employ gradients of the matrix elements with
respect to the optimization parameters.

These optimization techniques involve numerous evaluations of the Hamiltonian
matrix elements and their gradients. Therefore the analytica matrix elements
are of particular importance for the method.

##  3 Matrix elements

###  3.1 Overlap

The overlap $\langle g^{\prime}|g\rangle$ between a shifted Gaussian
$|g\rangle$ with parameters $A$, $\mathbf{s}$ and a shifted Gaussian
$|g^{\prime}\rangle$ with parameters $A^{\prime}$, $\mathbf{s}^{\prime}$ is
given as111 The overlap can be evaluated by an orthogonal coordinate
transformation, $\mathbf{r}=Q\mathbf{x}$, to the basis where the matrix $B$ is
diagonal: $B=QDQ^{\mathsf{T}}$ where $Q^{\mathsf{T}}Q=QQ^{\mathsf{T}}=1$ and
$D$ is a diagonal matrix,  $\displaystyle\langle g^{\prime}|g\rangle$
$\displaystyle=$ $\displaystyle\int d^{3}\vec{r}_{1}\dots
d^{3}\vec{r}_{N}\exp\left(-{\mathbf{r}}^{\mathsf{T}}B{\mathbf{r}}+\mathbf{v}^{\mathsf{T}}\mathbf{r}\right)=\int
d^{3}\vec{x}_{1}\dots
d^{3}\vec{x}_{N}\exp\left(-\sum_{i=1}^{N}{\vec{x}}_{i}\cdot
D_{ii}{\vec{x}}_{i}+\sum_{i=1}^{N}\vec{v}_{i}\cdot\vec{x}_{i}\right)$
$\displaystyle=$ $\displaystyle\prod_{i=1}^{N}\int
d^{3}\vec{x}_{i}\exp\left(-{\vec{x}}_{i}\cdot
D_{ii}{\vec{x}}_{i}+\vec{v}_{i}\cdot\vec{x}_{i}\right)=\prod_{i=1}^{N}\exp\left(\frac{1}{4D_{ii}}\vec{v}_{i}^{2}\right)\left(\frac{\pi}{D_{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}\;.$

$\langle
g^{\prime}|g\rangle=e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B)}\right)^{3/2}\equiv
M\,,\;$ (8)

where $B=A^{\prime}+A$, $\mathbf{v}=\mathbf{s}^{\prime}+\mathbf{s}$.

###  3.2 Kinetic energy

The non-relativistic kinetic energy operator $\hat{K}$ for an $N$-body system
of particles with coordinates $\vec{r}_{i}$ and masses $m_{i}$ is given as

$\hat{K}=-\sum_{i=1}^{N}\frac{\hbar^{2}}{2m_{i}}\frac{\partial^{2}}{\partial\vec{r}_{i}^{\,2}}\,.$
(9)

For completeness we shall 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}}}\,,$ (10)

where $\Lambda$ is a symmetric positive-definite matrix. The matrix element of
this operator is given as222 We first calculate two integrals,  $\left\langle
g^{\prime}\left|\mathbf{r}\right|g\right\rangle=\left(\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}=\mathbf{u}e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B)}\right)^{3/2}\,,$
(11) where $\mathbf{u}=\frac{1}{2}B^{-1}\mathbf{v}$, and
$\displaystyle\left\langle
g^{\prime}\left|\mathbf{r}^{\mathsf{T}}F\mathbf{r}\right|g\right\rangle=\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{\pi^{N}}{\det(B)}\right)^{3/2}=\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}\,,$
(12) from which the sought integral,  $\left\langle
g^{\prime}\left|-\frac{\partial}{\partial\bf
r}\Lambda\frac{\partial}{\partial\mathbf{r}^{\mathsf{T}}}\right|g\right\rangle=\left\langle
g^{\prime}\left|(\mathbf{s}^{\prime}-2A^{\prime}\mathbf{r})^{\mathsf{T}}\Lambda(\mathbf{s}-2A\mathbf{r})\right|g\right\rangle\,,$
(13) follows directly.

$\displaystyle\left\langle g^{\prime}\left|-\frac{\partial}{\partial\bf
r}\Lambda\frac{\partial}{\partial\mathbf{r}^{\mathsf{T}}}\right|g\right\rangle=\left\langle
g^{\prime}\left|(\mathbf{s}^{\prime\mathsf{T}}-2\mathbf{r}^{\mathsf{T}}A^{\prime})\Lambda(\mathbf{s}-2A\mathbf{r})\right|g\right\rangle$
$\displaystyle=\left(6\,\mathrm{trace}(A^{\prime}\Lambda
AB^{-1})+(\mathbf{s}^{\prime}-2A^{\prime}\mathbf{u})^{\mathsf{T}}\Lambda(\mathbf{s}-2A\mathbf{u})\right)M\,,$
(14)

where $\mathbf{u}\doteq\frac{1}{2}B^{-1}\mathbf{v}$, and $M$ is the overlap.

###  3.3 Potential energy

####  3.3.1 Central potential

A one-body central potential, $V(\vec{r}_{i})$, and a two-body central
potential, $V(\vec{r}_{i}-\vec{r}_{j})$, can be written in a convenient
general form, $V(w^{\mathsf{T}}\mathbf{r})$, where $w$ is a size-$N$ column of
numbers with all components equal zero except for $w_{i}=1$ for the one-body
potential and $w_{i}=-w_{j}=1$ for the two body potential.

##### Gaussian form-factor

For the Gaussian form-factor, $V(w^{\mathsf{T}}\mathbf{r})\propto
e^{-\gamma\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}}$, the matrix
element directly follows from the overlap integral,

$\left\langle
g^{\prime}\left|e^{-\gamma\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}}\right|g\right\rangle=e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{\prime-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B^{\prime})}\right)^{3/2}\equiv
M^{\prime}\,,$ (15)

where the matrix $B^{\prime}$ is a rank-1 update of the matrix
$B=A^{\prime}+A$, $B^{\prime}=B+\gamma ww^{\mathsf{T}}$.

If the determinant and the inverse of the matrix $B$ are known, their rank-1
updates can be calculated efficiently using the update formulas,

$\det(B+ab^{\mathsf{T}})=(1+b^{\mathsf{T}}B^{-1}a)\det(B)\,,$ (16)
$(B+ab^{\mathsf{T}})^{-1}=B^{-1}-{B^{-1}ab^{T}B^{-1}\over 1+b^{T}B^{-1}a}\,,$
(17)

where $a$ and $b$ are size-$N$ columns of numbers.

##### Oscillator form-factor

Another potential with a simple analytic matrix element is the oscillator
potential,
$V(w^{\mathsf{T}}\mathbf{r})\propto\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}$,
relevant for cold atoms in traps. The matrix element directly follows from
(12),

$\left\langle
g^{\prime}\left|\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}\right|g\right\rangle=\left(\frac{3}{2}w^{\mathsf{T}}B^{-1}w+\mathbf{u}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{u}\right)M\,.$
(18)

##### Other analytic form-factors

For a potential with a general form-factor, $V\propto
f(w^{\mathsf{T}}\mathbf{r})$, the matrix element reduces to a three-
dimensional integral,333 Suppose the form-factor $f(\vec{r})$ has a Fourier-
transform $\mathcal{F}(\vec{k})$, then  $\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\mathcal{F}(\vec{k})\left\langle
g^{\prime}\left|e^{i\vec{k}w^{\mathsf{T}}\mathbf{r}}\right|g\right\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}}\mathcal{F}(\vec{k})e^{-\alpha
k^{2}+i\vec{k}\vec{q}}$ (19) where $\alpha=\frac{1}{4}w^{\mathsf{T}}B^{-1}w$,
$\vec{q}=\frac{1}{2}w^{\mathsf{T}}B^{-1}\mathbf{v}$. Now the last integral can
as well be written as  $\left(\frac{\beta}{\pi}\right)^{\frac{3}{2}}\int
d^{3}\vec{r}\,f(\vec{r})\,e^{-\beta(\vec{r}-\vec{q})^{2}}\,,$ (20) where
$\beta=\frac{1}{4\alpha}=(w^{\mathsf{T}}B^{-1}w)^{-1}$.

$\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=M\left(\frac{\beta}{\pi}\right)^{\frac{3}{2}}\int
d^{3}\vec{r}\,f(\vec{r})\,e^{-\beta(\vec{r}-\vec{q})^{2}}\,,$ (21)

where $\beta=\left(w^{\mathsf{T}}B^{-1}w\right)^{-1}$ and
$\vec{q}=w^{\mathsf{T}}\mathbf{u}$.

If the potential does not depend on the direction of its argument the integral
reduces further to a one-dimensional integral,

$\displaystyle\left\langle
g^{\prime}\left|f(|w^{\mathsf{T}}\mathbf{r}|)\right|g\right\rangle$
$\displaystyle=$ $\displaystyle
M\left(\frac{\beta}{\pi}\right)^{\frac{3}{2}}2\pi\frac{e^{-\beta q^{2}}}{\beta
q}\int_{0}^{\infty}rdr\,f(r)\,e^{-\beta r^{2}}\sinh(2\beta qr)$
$\displaystyle\equiv$ $\displaystyle MJ[f]\,,$

where

$J[f]\doteq\left(\frac{\beta}{\pi}\right)^{\frac{3}{2}}2\pi\frac{e^{-\beta
q^{2}}}{\beta q}\int_{0}^{\infty}rdr\,f(r)\,e^{-\beta r^{2}}\sinh(2\beta
qr)\,.$ (23)

The integral (23) gives relatively simple analytic results for the Coulomb
form-factor, $1/r$,

$J\left[\frac{1}{r}\right]=\frac{\mathrm{erf}(\sqrt{\beta}q)}{q}\xrightarrow[q\to
0]{}\frac{2}{\sqrt{\pi}}\sqrt{\beta}\,,$ (24)

and for several short-range form-factors:

  * •

screened Yukawa, $e^{-\gamma r^{2}-\mu r}/r$,

$\displaystyle J\left[\frac{e^{-\gamma r^{2}-\mu r}}{r}\right]=$ (25)
$\displaystyle\frac{\sqrt{\beta}e^{-\beta
q^{2}}}{2q\sqrt{\beta+\gamma}}\left(e^{\frac{(\mu-2\beta
q)^{2}}{4(\beta+\gamma)}}\left(\mathrm{erf}\left(\frac{2\beta
q-\mu}{2\sqrt{\beta+\gamma}}\right)+1\right)+e^{\frac{(\mu+2\beta
q)^{2}}{4(\beta+\gamma)}}\left(\mathrm{erf}\left(\frac{\mu+2\beta
q}{2\sqrt{\beta+\gamma}}\right)-1\right)\right)\,,$

  * •

Yukawa, $e^{-\mu r}/r$,

$\displaystyle J\left[\frac{e^{-\mu r}}{r}\right]=\lim_{\gamma\to
0}J\left[\frac{e^{-\gamma r^{2}-\mu r}}{r}\right]=$ (26)
$\displaystyle\frac{e^{\frac{\mu^{2}}{4\beta}-\mu q}}{2q}\left(1-e^{2\mu
q}+\mathrm{erf}\left(\frac{2\beta q-\mu}{2\sqrt{\beta}}\right)+e^{2\mu
q}\mathrm{erf}\left(\frac{2\beta q+\mu}{2\sqrt{\beta}}\right)\right)\,,$ (27)

  * •

screened Coulomb, $e^{-\gamma r^{2}}/r$,

$\displaystyle J\left[\frac{e^{-\gamma r^{2}}}{r}\right]=\lim_{\mu\to
0}J\left[\frac{e^{-\gamma r^{2}-\mu
r}}{r}\right]=\frac{\sqrt{\beta}e^{-\frac{\beta\gamma
q^{2}}{\beta+\gamma}}\mathrm{erf}\left(\frac{\beta
q}{\sqrt{\beta+\gamma}}\right)}{q\sqrt{\beta+\gamma}}\,,$ (28)

  * •

screened exponential, $e^{-\gamma r^{2}-\mu r}$,

$\displaystyle J\left[e^{-\gamma r^{2}-\mu
r}\right]=\frac{-\partial}{\partial\mu}J\left[\frac{e^{-\gamma r^{2}-\mu
r}}{r}\right]=\frac{-\sqrt{\beta}}{4q(\beta+\gamma)^{3/2}}e^{\frac{(\mu-2\beta
q)^{2}}{4(\beta+\gamma)}-\beta q^{2}}$ (29) $\displaystyle\left((\mu-2\beta
q)\mathrm{erf}\left(\frac{2\beta
q-\mu}{2\sqrt{\beta+\gamma}}\right)+(\mu+2\beta q)e^{\frac{2\beta\mu
q}{\beta+\gamma}}\mathrm{erf}\left(\frac{\mu+2\beta
q}{2\sqrt{\beta+\gamma}}\right)+\mu-2\beta qe^{\frac{2\beta\mu
q}{\beta+\gamma}}-\mu e^{\frac{2\beta\mu q}{\beta+\gamma}}-2\beta q\right)\,,$

  * •

exponential, $e^{-\mu r}$,

$\displaystyle J\left[e^{-\mu r}\right]=\lim_{\gamma\to 0}J\left[e^{-\gamma
r^{2}-\mu r}\right]=\frac{-e^{\frac{\mu(\mu-4\beta q)}{4\beta}}}{4\beta q}$
(30) $\displaystyle\left((\mu-2\beta q)\mathrm{erf}\left(\frac{2\beta
q-\mu}{2\sqrt{\beta}}\right)+e^{2\mu q}(\mu+2\beta
q)\mathrm{erf}\left(\frac{\mu+2\beta q}{2\sqrt{\beta}}\right)+\mu-2\beta
qe^{2\mu q}-2\beta q-\mu e^{2\mu q}\right)\,.$

  * •

Gaussian, $e^{-\gamma r^{2}}$,

$J\left[e^{-\gamma r^{2}}\right]=\lim_{\mu\to 0}J\left[e^{-\gamma r^{2}-\mu
r}\right]=\left(\frac{\beta}{\beta+\gamma}\right)^{3/2}e^{-\frac{\beta\gamma
q^{2}}{\beta+\gamma}}\,.$ (31)

Since the error function has analytic derivative,
$\mathrm{erf}^{\prime}(x)=2e^{-x^{2}}/\sqrt{\pi}$, the following form-factors
also have analytic matrix elements,

$J\left[r^{n}\frac{e^{-\gamma r^{2}-\mu
r}}{r}\right]=\left(\frac{-\partial}{\partial\mu}\right)^{n}J\left[\frac{e^{-\gamma
r^{2}-\mu r}}{r}\right]\,.$ (32)

Although there are few other form-factors with analytic matrix elements, the
resulting relatively complicated expressions are not conducive to subsequent
analytic calculations of gradients with respect to the elements of matrices
$A$ and the shift-vectors $\mathbf{s}$ through the quantities $\beta$ and $q$.
In such cases it is probably more efficient to represent the potential as a
linear combinations of several Gaussians.

####  3.3.2 Tensor potential

The tensor potential $V_{\mathrm{t}}$ between two particles with coordinates
$\vec{r}_{1}$, $\vec{r}_{2}$ and spins $\vec{S}_{1}$, $\vec{S}_{2}$ can be
written in the form [11]

$V_{\mathrm{t}}(r)\propto
f(r)(\vec{S}_{1}\cdot\vec{r})(\vec{S}_{2}\cdot\vec{r})\,,$ (33)

where $\vec{r}=\vec{r}_{1}-\vec{r}_{2}$ is the relative coordinate between the
particles, and $f(r)$ is the radial form-factor of the potential. In this form
the potential has a central spin-spin component,

$\frac{1}{3}f(r)r^{2}(\vec{S}_{1}\cdot\vec{S}_{2})\,,$ (34)

which is often subtracted from the above form to make sure the potential
contains only the spherical tensor component.

Introducing the size-$N$ column of numbers $w$,

$w=\\{w_{i}|w_{1}=1,w_{2}=-1,w_{i\neq 1,2}=0\\}\,,$ (35)

and vector-columns $\mathbf{y}_{1}=\vec{S}_{1}w$,
$\mathbf{y}_{2}=\vec{S}_{2}w$, the tensor potential can be written in a
convenient general form,

$\hat{V}_{\mathrm{t}}\propto
f(w^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{r})\,,$
(36)

where

$\mathbf{y}^{\mathsf{T}}\mathbf{r}\equiv\sum_{i=1}^{N}\vec{y}_{i}\cdot\vec{r}_{i}\,.$
(37)

The tensor matrix element can be represented as a derivative of the central
matrix element with the same form-factor,

$\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=\left(\mathbf{y}_{1}^{\mathsf{T}}\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\left(\mathbf{y}_{2}^{\mathsf{T}}\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle\,.$ (38)

The central part (36) of the tensor potential has a similar analytic
representation,

$\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})(w^{\mathsf{T}}\mathbf{r})(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=\left(w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\left(w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle\,.$ (39)

Thus if the matrix element of the central potential with a given form-factor
$f(w^{\mathsf{T}}\mathbf{r})$ is analytic, the matrix element of the tensor
potential
$f(w^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{r})$
is also analytic. In particular, the tensor matrix elements of the screened
Yukawa form-factor and its descendants (32) are also analytic.

##### Gaussian form-factor

For a Gaussian form-factor the matrix elements (38,39) are readily given as

$\displaystyle\left\langle
g^{\prime}\left|e^{-\gamma\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}}(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=\left(\frac{1}{2}\mathbf{y}_{1}^{\mathsf{T}}B^{\prime-1}\mathbf{y}_{2}+(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{u}^{\prime})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{u}^{\prime})\right)M^{\prime}\,,$
$\displaystyle\left\langle
g^{\prime}\left|e^{-\gamma\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r}}(\mathbf{r}^{\mathsf{T}}ww^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=\left(\frac{3}{2}w^{\mathsf{T}}B^{\prime-1}w+\mathbf{u}^{\prime\mathsf{T}}ww^{\mathsf{T}}\mathbf{u}^{\prime}\right)M^{\prime}\;,$

where $B^{\prime}=B+\gamma ww^{\mathsf{T}}$, $B=A^{\prime}+A$,
$\mathbf{u}^{\prime}=\frac{1}{2}B^{\prime-1}\mathbf{v}$,
$\mathbf{v}=\mathbf{s}^{\prime}+\mathbf{s}$, and where

$M^{\prime}=e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{\prime-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B^{\prime})}\right)^{3/2}$
(40)

is the central Gaussian matrix element.

Again the updates $\det(B+\gamma ww^{\mathsf{T}})$ and $(B+\gamma
ww^{\mathsf{T}})^{-1}$ can be efficiently calculated using rank-1 update
formulas.

##### Other form-factors

The tensor matrix element (38) can be written in component form as

$\displaystyle\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{1}^{\mathsf{T}}\mathbf{r})(\mathbf{y}_{2}^{\mathsf{T}}\mathbf{r})\right|g\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{i,j=1}^{N}\sum_{a,b=1}^{3}(\vec{y}_{1})_{ia}(\vec{y}_{2})_{jb}\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\vec{r}_{ia}\vec{r}_{jb}\right|g\right\rangle$
(41) $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{N}\sum_{a,b=1}^{3}(\vec{y}_{1})_{ia}(\vec{y}_{2})_{jb}\frac{\partial}{\partial\vec{v}_{ia}}\frac{\partial}{\partial\vec{v}_{jb}}\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle\,.$

where $\vec{r}_{ia}$ is the number-$a$ component of the vector $\vec{r}_{i}$,
and where

$\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=MJ\,,$ (42)
$M=e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B)}\right)^{\frac{3}{2}}\,,$
(43) $J=\left(\frac{\beta}{\pi}\right)^{\frac{3}{2}}2\pi\frac{e^{-\beta
q^{2}}}{\beta q}\int_{0}^{\infty}rdr\,f(r)\,e^{-\beta r^{2}}\sinh(2\beta
qr)\,.$ (44)

The derivative of the central matrix element at the right-hand side of (41)
can now be evaluated as

$\displaystyle\frac{\partial}{\partial\vec{v}_{ia}}\frac{\partial}{\partial\vec{v}_{jb}}\left(MJ\right)=\left(\frac{\partial}{\partial\vec{v}_{ia}}\frac{\partial}{\partial\vec{v}_{jb}}M\right)J+M\left(\frac{\partial}{\partial\vec{v}_{ia}}\frac{\partial}{\partial\vec{v}_{jb}}J\right)$
$\displaystyle+\left(\frac{\partial}{\partial\vec{v}_{ia}}M\right)\left(\frac{\partial}{\partial\vec{v}_{jb}}J\right)+\left(\frac{\partial}{\partial\vec{v}_{jb}}M\right)\left(\frac{\partial}{\partial\vec{v}_{ia}}J\right)\,,$
(45)

where

$\frac{\partial}{\partial\vec{v}_{ia}}M=\sum_{k=1}^{N}\frac{1}{2}B^{-1}_{ik}\vec{v}_{ka}M\,,$
(46)
$\frac{\partial}{\partial\vec{v}_{jb}}\frac{\partial}{\partial\vec{v}_{ia}}M=\frac{1}{2}B^{-1}_{ji}\delta_{ab}M+\sum_{k,l=1}^{N}\frac{1}{2}B^{-1}_{ik}\vec{v}_{ka}\frac{1}{2}B^{-1}_{jl}\vec{v}_{lb}M\,,$
(47) $\frac{\partial}{\partial\vec{v}_{ia}}J=\frac{\partial J}{\partial
q}\frac{\partial q}{\partial\vec{v}_{ia}}\,,$ (48)
$\frac{\partial}{\partial\vec{v}_{jb}}\frac{\partial}{\partial\vec{v}_{ia}}J=\frac{\partial^{2}J}{\partial
q^{2}}\frac{\partial q}{\partial\vec{v}_{jb}}\frac{\partial
q}{\partial\vec{v}_{ia}}+\frac{\partial J}{\partial
q}\frac{\partial^{2}q}{\partial\vec{v}_{jb}\partial\vec{v}_{ia}}\,,$ (49)
$\frac{\partial
q}{\partial\vec{v}_{ia}}=\frac{1}{q}\sum_{k=1}^{N}h_{i}h_{k}\vec{v}_{ka}\,,$
(50)
$\frac{\partial^{2}q}{\partial\vec{v}_{jb}\partial\vec{v}_{ia}}=\frac{1}{q}h_{i}h_{j}\delta_{ab}-\frac{2}{q^{2}}\sum_{k=1}^{N}h_{i}h_{k}\vec{v}_{ka}\sum_{l=1}^{N}h_{j}h_{l}\vec{v}_{lb}\,.$
(51)

where $h=w^{\mathsf{T}}\frac{1}{2}B^{-1}$.

Now to finish the calculation one only needs to calculate $\partial J/\partial
q$ and $\partial^{2}J/\partial q^{2}$ for the given potential. For the
screened Yukawa form-factor and its descendants (32) the derivatives $\partial
J/\partial q$ and $\partial^{2}J/\partial q^{2}$ are analytic although the
actual calculations are relatively tedious and should be best performed by a
computer algebra software like Maxima [12]. For example, the expression for

$\frac{\partial^{2}}{\partial q^{2}}J\left[\frac{e^{-\gamma r^{2}-\mu
r}}{r}\right]$ (52)

can be readily obtained by the following Maxima script,

    assume(beta>0,gamma>0,mu>0,q>0);
    J(f) := (beta/%pi)^(3/2)*2*%pi/beta/q*exp(-beta*q^2)
    *integrate(r*f*exp(-beta*r^2)*sinh(2*beta*q*r),r,0,inf);
    fortran(diff(J(exp(-gamma*r^2-mu*r)/r), q, 2));
    Ψ

which analytically calculates $\frac{\partial^{2}J}{\partial q^{2}}$ for the
screened Yukawa potential and outputs the corresponding Fortran code.

Even if $J$ is not analytic one can evaluate its derivatives by calculating
numerically a few extra integrals,

$\displaystyle\frac{\partial J}{\partial
q}=-{{2\,\sqrt{\beta}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q^{2}}}$
$\displaystyle-{{4\,\beta^{{{3}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}}}$
$\displaystyle+{{4\,\beta^{{{3}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r^{2}\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\cosh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q}}\,,$
(53) $\displaystyle\frac{\partial^{2}J}{\partial
q^{2}}={{8\,\beta^{{{5}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r^{3}\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q}}$
$\displaystyle+{{8\,\beta^{{{5}\over{2}}}\,q\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}}}$
$\displaystyle+{{4\,\beta^{{{3}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q}}$
$\displaystyle+{{4\,\sqrt{\beta}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\sinh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q^{3}}}$
$\displaystyle-{{8\,\beta^{{{3}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r^{2}\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\cosh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}\,q^{2}}}$
$\displaystyle-{{16\,\beta^{{{5}\over{2}}}\,e^{-\beta\,q^{2}}\,\int_{0}^{\infty}{r^{2}\,e^{-\beta\,r^{2}}\,f\left(r\right)\,\cosh\left(2\,\beta\,q\,r\right)\;dr}}\over{\sqrt{\pi}}}\,.$
(54)

####  3.3.3 Spin-orbit potential

The spin-orbit potential between two particles with coordinates $\vec{r}_{1}$
and $\vec{r}_{2}$ and spins $\vec{S}_{1}$ and $\vec{S}_{2}$ can be written in
the form [11]

$V_{\mathsf{so}}\propto f(r)\left(\vec{S}\cdot\vec{L}\right)\,,$ (55)

where $\vec{r}=\vec{r}_{1}-\vec{r}_{2}$ is the relative coordinate between the
particles; $f(r)$ is the radial form-factor of the potential;
$\vec{S}=\vec{S}_{1}+\vec{S}_{2}$ is the total spin of the two particles; and
$\vec{L}$ is the relative orbital momentum between the two particles,

$\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)\,,$
(56)

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

The orbital momentum operator can be written, using the size-$N$ column of
numbers $w$,

$w=\\{w_{i}|w_{1}=1,w_{2}=-1,w_{i\neq 1,2}=0\\}\,,$ (57)

in the 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)\,,$ (58)

where

$w^{\mathsf{T}}\mathbf{r}\equiv\sum_{i=1}^{N}w_{i}\vec{r}_{i}\,,\;w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{r}^{\mathsf{T}}}\equiv\sum_{i=1}^{N}w_{i}\frac{\partial}{\partial\vec{r}_{i}}\,.$
(59)

For a given form-factor $f(w^{\mathsf{T}}\mathbf{r})$ the spin-orbit matrix
element can be represented through the central matrix element with the same
form-factor,

$\displaystyle\left\langle
g^{\prime}\left|f(w^{\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=\left(w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\times
w^{\mathsf{T}}\left(\mathbf{s}-2A\frac{\partial}{\partial\mathbf{v}^{\mathsf{T}}}\right)\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle\,,$ (60)

so that if the central matrix element is analytic — as is the case for the
screened Yukawa form-factor and its descendants (32) — the spin-orbit matrix
element is also analytic.

##### Gaussian form-factor

For a Gaussian form-factor the spin-orbit matrix element is given as

$\displaystyle\left\langle
g^{\prime}\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=\frac{1}{2}w^{\mathsf{T}}B^{\prime-1}\mathbf{v}\times\left(w^{\mathsf{T}}\mathbf{s}-w^{\mathsf{T}}AB^{\prime-1}\mathbf{v}\right)M^{\prime}\,,$
(61)

where $B^{\prime}=B+\gamma ww^{\mathsf{T}}$, $B=A^{\prime}+A$,
$\mathbf{u}^{\prime}=\frac{1}{2}B^{\prime-1}\mathbf{v}$,
$\mathbf{v}=\mathbf{s}^{\prime}+\mathbf{s}$, and where

$M^{\prime}=e^{\frac{1}{4}\mathbf{v}^{\mathsf{T}}B^{\prime-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B^{\prime})}\right)^{3/2}$
(62)

is the central Gaussian matrix element.

Again $\det(B^{\prime})$ and ${B^{\prime}}^{-1}$ can be efficiently calculated
using rank-1 update formulas.

##### General form-factor

The component number $a$ of the spin-orbit matrix element (60) can be written
in component form as

$\displaystyle\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\left(w^{\mathsf{T}}\mathbf{r}\times
w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{r}^{\mathsf{T}}}\right)_{a}\right|g\right\rangle=$
(63)
$\displaystyle\sum_{b,c=1}^{3}\epsilon_{abc}\sum_{k=1}^{N}w_{k}\frac{\partial}{\partial\vec{v}_{kb}}\sum_{l=1}^{N}w_{l}\left(\vec{s}_{lc}-2\sum_{j=1}^{N}A_{lj}\frac{\partial}{\partial\vec{v}_{jc}}\right)\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle\,.$ (64)

It is clearly a linear combination of the first,
$\frac{\partial}{\partial\vec{v}_{kb}}\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle$, and second,
$\frac{\partial}{\partial\vec{v}_{kb}}\frac{\partial}{\partial\vec{v}_{jc}}\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle$, derivatives
of the corresponding central matrix element $\left\langle
g^{\prime}\left|f(w^{\mathsf{T}}\mathbf{r})\right|g\right\rangle=MJ$. These
quantities have been calculated in the previous chapter.

###  3.4 Many-body forces

Many-body potentials in nuclear physics have form-factors which depend on the
coordinates of several nucleons. In addition they may have tensor and spin-
orbit post-factors.

For a many-body potential with a Gaussian form-factor,
$e^{-\mathbf{r}^{\mathsf{T}}W\mathbf{r}}$, where $W$ is a symmetric positive-
definite matrix, the central matrix element is given as

$\left\langle
g^{\prime}\left|e^{-\mathbf{r}^{\mathsf{T}}W\mathbf{r}}\right|g\right\rangle=e^{\frac{1}{4}\mathbf{v}(B+W)^{-1}\mathbf{v}}\left(\frac{\pi^{N}}{\det(B+W)}\right)^{3/2}\,.$
(65)

The simple tensor and spin-orbit post-factors can be calculated in the same
way as has been done in the previous two chapters.

##  4 Conclusion

In quantum few-body physics the analytic matrix elements are of importance for
the correlated Gaussians method as they facilitate extensive numerical
optimizations of the variational wave-functions. In this paper it has been
shown that potentials in the form

$f(w^{\mathsf{T}}\mathbf{r})\,,\;f(w^{\mathsf{T}}\mathbf{r})(\mathbf{y}^{\mathsf{T}}\mathbf{r})(\mathbf{z}^{\mathsf{T}}\mathbf{r})\,,\;f(w^{\mathsf{T}}\mathbf{r})(a^{\mathsf{T}}\mathbf{r})(b^{\mathsf{T}}\mathbf{r})\,,\;f(w^{\mathsf{T}}\mathbf{r})\left(w^{\mathsf{T}}\mathbf{r}\times
w^{\mathsf{T}}\frac{\partial}{\partial\mathbf{r}^{\mathsf{T}}}\right)\,,$ (66)

have analytic matrix elements between shifted correlated Gaussians for the
following class of form-factors,

$f(r)=r^{n}\frac{e^{-\gamma r^{2}-\mu r}}{r}\,,\;n=0,1,\dots\,.$ (67)

Of these form-factors the Gaussian form-factors produce particularly simple
and concise analytic expressions. Therefore an efficient strategy could be to
represent the potentials at hand as linear combinations of Gaussians and then
use the analytic expressions for the Gaussians.

## References

  * [1] J. Mitroy et al., Theory and applications of explicitly correlated Gaussians, Rev. Mod. Phys. 85, 693 (2013).
  * [2] Sergiy Bubin, Michele Pavanello, Wei-Cheng Tung, Keeper L. Sharkey, and Ludwik Adamowicz, Born-Oppenheimer and Non-Born-Oppenheimer Atomic and Molecular Calculations with Explicitly Correlated Gaussians, Chemical Reviews 113, 36 (2013).
  * [3] Mauricio Cafiero, Ludwik Adamowicz, Analytical Gradients for Singer’s Multicenter n-Electron Explicitly Correlated Gaussians, International Journal of Quantum Chemistry, Vol. 82, 151–159 (2001).
  * [4] Keeper L. Sharkey, Sergiy Bubin, and Ludwik Adamowicz, An algorithm for calculating atomic D states with explicitly correlated Gaussian functions, The Journal of Chemical Physics 134, 044120 (2011).
  * [5] K. M. Daily and Chris H. Greene, Extension of the correlated Gaussian hyperspherical method to more particles and dimensions, Physical Review A 89, 012503 (2014)
  * [6] Wei-Cheng Tung and Ludwik Adamowicz, Accurate potential energy curve of the LiH+ molecule calculated with explicitly correlated Gaussian functions, J. Chem. Phys. 140, 124315 (2014).
  * [7] Sergiy Bubin, Martin Formanek, Ludwik Adamowicz, Universal all-particle explicitly-correlated Gaussians for non-Born–Oppenheimer calculations of molecular rotationless states, Chemical Physics Letters 647, 122 (2016).
  * [8] X.Y. Yin and D. Blume, Trapped unitary two-component Fermi gases with up to ten particles, Physical Review A 92, 013608 (2015).
  * [9] M. Galassi et al, GNU Scientific Library Reference Manual - Third Edition (January 2009), ISBN 0954612078.
  * [10] Y. Suzuki and K. Varga, Stochastic Variational Approach to Quantum-Mechanical Few-Body Problems, ISBN 3-540-65152-7, Springer-Verlag, Berlin, 1998.
  * [11] Peter Ring and Peter Schuck, The Nuclear Many-Body Problem, ISBN 3-540-09820-8, Springer-Verlag, Berlin, 1980.
  * [12] Maxima.sourceforge.net. Maxima, a Computer Algebra System. Version 5.34.1 (2014). http://maxima.sourceforge.net/