Path integral molecular dynamics approximations of quantum canonical observables

Xin Huang Institutionen för Matematik, Kungl. Tekniska Högskolan, 100 44 Stockholm, Sweden [email protected] , Petr Plecháč Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA [email protected] , Mattias Sandberg Institutionen för Matematik, Kungl. Tekniska Högskolan, 100 44 Stockholm, Sweden [email protected] and Anders Szepessy Institutionen för Matematik, Kungl. Tekniska Högskolan, 100 44 Stockholm, Sweden [email protected]

Abstract.

Mean-field molecular dynamics based on path integrals is used to approximate canonical quantum observables for particle systems consisting of nuclei and electrons. A computational bottleneck is the sampling from the Gibbs density of the electron operator, which due to the fermion sign problem has a computational complexity that scales exponentially with the number of electrons. In this work we construct an algorithm that approximates the mean-field Hamiltonian by path integrals for fermions. The algorithm is based on the determinant of a matrix with components based on Brownian bridges connecting permuted electron coordinates. The computational work for $n$ electrons is $\mathcal{O}(n^{3})$ , which reduces the computational complexity associated with the fermion sign problem. We analyze a bias resulting from this approximation and provide a computational error indicator. It remains to rigorously explain the surprisingly high accuracy.

1. Background to approximations of quantum canonical observables

Path integral methods are used to approximate observables in the canonical quantum ensemble of particle systems consisting of nuclei and electrons. Specific implementations include, for instance, ring polymer molecular dynamics, centroid molecular dynamics, and related mean-field molecular dynamics, see [27] and [20]. A computational challenge arises especially for systems of fermions from the so-called (fermion) sign problem, causing the computational work of Monte Carlo approximations to grow exponentially with the number of fermion particles for a fixed expected approximation error [11, 25]. The aim of this work is to provide an alternative path integral formulation that reduces the sign problem, where the Monte Carlo method averages the determinant of a matrix with components based on Brownian bridge path integrals. The formulation uses the Feynman-Kac representation of the Gibbs partition function for the electron operator, with permuted Brownian bridge end coordinates due to indistinguishable particles. In this setting the partition function is represented by a sum over permutations of a rank four tensor. Sums of permutations of tensors are in general computationally feasible only in low dimension, since such sums often are NP-hard, see [18], with the exception of computing determinants of matrices. To reduce the tensor problem to determinants we use that the Gibbs potential energy for a nuclei-electron system can be written as a product based on path integrals for the energy with respect to each particle. In other words, we use that the potential energy is based on Coulomb pair interactions and its symmetry with respect to permutations of particles. Our formulation has no bias in the case of a separable potential, i.e., when it can be written as a sum of potentials that each depends on one particle only. Numerical tests show that the approximation error is relatively small also for non-separable potentials based on Coulomb electron repulsion, since the marginal distribution of initial and end particle positions turns out to be Gaussians.

We study the specific setting in which the electron kinetic energy operator is the Laplacian, which provides the Brownian bridge paths. Hence we consider the Hamiltonian

\widehat{H}=-\frac{1}{2M}\Delta_{\mathbf{X}}-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X})=:-\frac{1}{2M}\Delta_{\mathbf{X}}+H_{e}

with the Coulomb potential $V:\mathbb{R}^{3n}\times\mathbb{R}^{3N}\to\mathbb{R}$ , the nuclei coordinates $\mathbf{X}=(\mathbf{X}_{1},\ldots,\mathbf{X}_{N})\in\mathbb{R}^{3N}$ , and the electron coordinates $\mathbf{x}=(\mathbf{x}_{1},\ldots\mathbf{x}_{n})\in\mathbb{R}^{3n}$ , where $M$ is the nuclei-electron mass ratio. We use the Hartree atomic units where the reduced Planck constant $\hbar=1$ , the electron charge $e=1$ , the Bohr radius $a_{0}=1$ and the electron mass $m_{e}=1$ . The small semiclassical parameter $1/M$ is, for example, in the case of a proton-electron system $1/M=m_{e}/m_{p}\approx 1/1836$ . The objective is to approximate canonical quantum correlation observables based on the Hamiltonian $\widehat{H}$ and the inverse temperature $\beta>0$

\frac{\mathrm{Tr}\,(\widehat{A}_{t}\widehat{B}_{0}e^{-\beta\widehat{H}})}{\mathrm{Tr}\,(e^{-\beta\widehat{H}})}

and its symmetrized form

\mathfrak{T}_{{\rm qm}}:=\frac{\mathrm{Tr}\,\big{(}(\widehat{A}_{t}\widehat{B}_{0}+\widehat{B}_{0}\widehat{A}_{t})e^{-\beta\widehat{H}}\big{)}}{2\mathrm{Tr}\,(e^{-\beta\widehat{H}})}

for quantum observables $\widehat{A}_{t}=e^{{\rm i}t\sqrt{M}\widehat{H}}\widehat{A}_{0}e^{-{\rm i}t\sqrt{M}\widehat{H}}$ and $\widehat{B}_{0}$ , at times $t$ and $0$ , using the trace, $\mathrm{Tr}\,$ , of a quantum operator over the nuclei and electron degrees of freedom. Such correlation observables are used, for instance, to determine the diffusion constant, reaction rates and other transport coefficients, [32].

Our more precise aim is to obtain computable classical molecular dynamics approximations of quantum observables. For this purpose we need the Weyl symbol $A:\mathbb{R}^{6N}\to\mathbb{C}$ related to the quantum operator $\widehat{A}$ , defined by

(1.1)

\begin{split}\widehat{A}\phi(\mathbf{X})=&\int_{\mathbb{R}^{3N}}(\frac{\sqrt{M}}{2\pi})^{3N}\int_{\mathbb{R}^{3N}}e^{{\rm i}M^{1/2}(\mathbf{X}-\mathbf{Y})\cdot\mathbf{P}}A\big{(}\frac{1}{2}(\mathbf{X}+\mathbf{Y}),\mathbf{P}\big{)}{\mathrm{d}}\mathbf{P}\,\phi(\mathbf{Y}){\mathrm{d}}\mathbf{Y}\\ \end{split}

for Schwartz functions $\phi:\mathbb{R}^{3N}\to\mathbb{C}$ , see [33]. We assume that the Weyl symbols $A_{0}:\mathbb{R}^{3N}\times\mathbb{R}^{3N}\to\mathbb{C}$ and $B_{0}:\mathbb{R}^{3N}\times\mathbb{R}^{3N}\to\mathbb{C}$ for the initial quantum observables are independent of the electron coordinates. We will use that the canonical quantum observables can be approximated by mean-field molecular dynamics

(1.2)

\mathfrak{T}_{\mathrm{md}}(t):=\frac{\int_{\mathbb{R}^{6N}}A_{0}\big{(}\mathbf{Z}_{t}(\mathbf{Z}_{0})\big{)}B_{0}(\mathbf{Z}_{0})\mathrm{Tr}\,(e^{-\beta H(\mathbf{Z}_{0})}){\mathrm{d}}\mathbf{Z}_{0}}{\int_{\mathbb{R}^{6N}}\mathrm{Tr}\,(e^{-\beta H(\mathbf{Z}_{0})}){\mathrm{d}}\mathbf{Z}_{0}}

where the Hamiltonian is given by

H=\frac{|\mathbf{P}|^{2}}{2}-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X})=\frac{|\mathbf{P}|^{2}}{2}+H_{e}(\mathbf{x},\mathbf{X})

and the trace, $\mathrm{Tr}\,$ , of a Weyl symbol is the trace with respect to the electron degrees of freedom, that is the trace on an appropriate antisymmetric subset of $L^{2}(\mathbb{R}^{3n})$ . The phase-space nuclei coordinates $\mathbf{Z}_{t}:=(\mathbf{X}_{t},\mathbf{P}_{t})$ solve the mean-field Hamiltonian system

(1.3)

\begin{split}\dot{\mathbf{X}}_{t}=&\;\;\;\;\nabla_{\mathbf{P}}h(\mathbf{X}_{t},\mathbf{P}_{t})\,,\\ \dot{\mathbf{P}}_{t}=&-\nabla_{\mathbf{X}}h(\mathbf{X}_{t},\mathbf{P}_{t})\,,\\ \end{split}

with the initial data $\mathbf{Z}_{0}:=(\mathbf{X}_{0},\mathbf{P}_{0})\in\mathbb{R}^{3N}\times\mathbb{R}^{3N}$ of nuclei positions and momenta, where the mean-field Hamiltonian $h:\mathbb{R}^{3N}\times\mathbb{R}^{3N}\to\mathbb{R}$ is defined by

(1.4)

h(\mathbf{Z}):=\frac{\mathrm{Tr}\,(H(\mathbf{Z})e^{-\beta H(\mathbf{Z})})}{\mathrm{Tr}\,(e^{-\beta H(\mathbf{Z})})}\,.

Assuming that the electron problem for a given nuclei position,

H_{e}(\cdot,\mathbf{X})\psi_{i}(\cdot,\mathbf{X})=\lambda_{i}(\mathbf{X})\psi_{i}(\cdot,\mathbf{X})\,,

has the eigenvalues $\lambda_{i}(\mathbf{X})\in\mathbb{R}$ , and eigenfunctions $\psi_{i}(\cdot,\mathbf{X})\in L^{2}(\mathbb{R}^{3n})$ , $i=1,2,3,\ldots$ , then

h(\mathbf{X},\mathbf{P})=\frac{|\mathbf{P}|^{2}}{2}+\frac{\sum_{i=1}^{\infty}\lambda_{i}(\mathbf{X})e^{-\beta\lambda_{i}(\mathbf{X})}}{\sum_{i=1}^{\infty}e^{-\beta\lambda_{i}(\mathbf{X})}}=:\frac{|\mathbf{P}|^{2}}{2}+\lambda_{*}(\mathbf{X})\,.

The work [20] derives the approximation error

(1.5)

|\mathfrak{T}_{\mathrm{qm}}(t)-\mathfrak{T}_{\mathrm{md}}(t)|=\mathcal{O}(M^{-1}+t\epsilon_{1}^{2}+t^{2}\epsilon_{2}^{2})\,,

where

\begin{split}\epsilon_{1}^{2}&=\frac{\|\mathrm{Tr}\,\big{(}(H-h)^{2}e^{-\beta H}\big{)}\|_{L^{1}(\mathbb{R}^{6N})}}{\|\mathrm{Tr}\,\big{(}e^{-\beta H}\big{)}\|_{L^{1}(\mathbb{R}^{6N})}}=\frac{\|\sum_{i=1}^{\infty}(\lambda_{i}-\lambda_{*})^{2}e^{-\beta\lambda_{i}}\|_{L^{1}(\mathbb{R}^{3N})}}{\|\sum_{i=1}^{\infty}e^{-\beta\lambda_{i}}\|_{L^{1}(\mathbb{R}^{3N})}}\,,\\ \epsilon_{2}^{2}&=\frac{\|\sum_{i=1}^{\infty}|\nabla(\lambda_{i}-\lambda_{*})|^{2}e^{-\beta\lambda_{i}}\|_{L^{1}(\mathbb{R}^{3N})}}{\|\sum_{i=1}^{\infty}e^{-\beta\lambda_{i}}\|_{L^{1}(\mathbb{R}^{3N})}}\,,\end{split}

for the case that the electron part, $-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X})$ , is approximated by a finite dimensional matrix. The mean-field Hamiltonian has the representation

(1.6)

\begin{split}h(\mathbf{X},\mathbf{P})&=\frac{|\mathbf{P}|^{2}}{2}+\frac{\mathrm{Tr}\,\Big{(}\big{(}-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X})\big{)}e^{-\beta(-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X}))}\Big{)}}{\mathrm{Tr}\,\big{(}e^{-\beta(-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X}))}\big{)}}\\ &=\frac{|\mathbf{P}|^{2}}{2}+\frac{\mathrm{Tr}\,(H_{e}e^{-\beta H_{e}})}{\mathrm{Tr}\,(e^{-\beta H_{e}})}=\frac{|\mathbf{P}|^{2}}{2}-\partial_{\beta}\log\big{(}\mathrm{Tr}\,(e^{-\beta H_{e}})\big{)}\,.\end{split}

To apply the molecular dynamics method (1.2) requires approximate sampling from the canonical Gibbs measure

\mathbf{X}\mapsto\mathrm{Tr}\,(e^{-\beta H_{e}(\cdot,\mathbf{X})})/\int_{\mathbb{R}^{3N}}\mathrm{Tr}\,(e^{-\beta H_{e}(\cdot,\mathbf{X}_{0})}){\mathrm{d}}\mathbf{X}_{0}:\mathbb{R}^{3N}\to\mathbb{R}\,.

In particular, approximation of the mean-field $h:\mathbb{R}^{6N}\to\mathbb{R}$ is needed. The purpose of this work is to formulate such approximations based on Monte Carlo sampled path integrals for fermion systems. The mean-field formulation (1.2) is closely related to the ring polymer molecular dynamics and centroid molecular dynamics methods, see [17, 27, 26].

The main result in this work is the construction and numerical tests of a proposed new computational method that approximates path integrals for fermions which is based on Brownian bridge paths and estimation of a suitable determinant in Section 3. We present the necessary background on path integrals in Section 2. In Section 3.1 we describe a perturbation analysis which, in combination with numerical experiments in Section 4, provides a computational error indicator to roughly estimate the accuracy of the developed determinant formulation. The numerical implementation and computational experiments for test cases with varying number of particles and inverse temperature are reported in Section 4. Finally, in Section 5 we comment on a surprisingly high accuracy for which it remains to find a rigorous theoretical answer.

2. Path integrals

First we present a brief background on path integral methods for the mean-field (1.6), to be used in Section 3 for the new representation of the partition function based on a determinant and Brownian bridge coordinate paths.

We denote $\mathbf{W}:[0,\beta]\to\mathbb{R}^{dn}$ the standard Wiener process with independent components and define a Wiener path starting at $\mathbf{x}_{0}$ to be $\mathbf{x}_{t}=\mathbf{x}_{0}+\mathbf{W}_{t}$ .

2.1. Distinguishable particles

If electrons are treated as distinguishable particles the Feynman-Kac formulation of path integrals, see [21, Theorem 7.6] yields the expected value representation

(2.1)

e^{-\beta H_{e}}\phi(\mathbf{x}_{0})=\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\phi(\mathbf{x}_{\beta})]\,.

Then the $L^{2}(\mathbb{R}^{dn})$ trace of the Gibbs density operator becomes

(2.2)

\mathrm{Tr}\,(e^{-\beta H_{e}(\cdot,\mathbf{X})})=\int_{\mathbb{R}^{dn}}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\delta(\mathbf{x}_{\beta}-\mathbf{x}_{0})]{\mathrm{d}}\mathbf{x}_{0}\,.

To approximate the mean-field Hamiltonian we compute

\begin{split}\frac{\mathrm{Tr}\,(H_{e}e^{-\beta H_{e}})}{\mathrm{Tr}\,(e^{-\beta H_{e}})}&=-\partial_{\beta}\log\mathrm{Tr}\,(e^{-\beta H_{e}})\\ &=-\partial_{\beta}\log\big{(}\frac{1}{(2\pi\beta)^{3n/2}}\int_{\mathbb{R}^{3n}}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\,|\,\mathbf{x}_{\beta}=\mathbf{x}_{0}]{\mathrm{d}}\mathbf{x}_{0}\big{)}\end{split}

where the conditional expectation is on Brownian bridge paths $\mathbf{x}_{t}$ , $\mathbf{x}_{\beta}=\mathbf{x}_{0}$ . Hence, the integral $\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t$ and the Brownian bridge $\mathbf{x}_{t}$ need to be discretized and generated. One alternative is to directly generate discrete approximations of Brownian bridges, which is the focus in this work. Another more common method, cf. [9, 11], is to discretize the path integral as follows. Let $\mathbf{x}(j):=\mathbf{x}(j\beta/J)$ be the values of the paths at the time steps $t_{j}:=j\beta/J=:j\Delta t$ . The normal distribution of the Brownian increments in the left hand side of (2.2) implies that the Euler approximation representation, for the case with distinguishable particles

(2.3)

\begin{split}&\mathrm{Tr}\,(e^{-\beta H_{e}(\cdot,\mathbf{X})})\\ &=\lim_{\Delta t\to 0+}\int_{\mathbb{R}^{dn}}\mathbb{E}[e^{-\sum_{j=0}^{J-1}V(\mathbf{x}(j),\mathbf{X}){\Delta t}}\delta\big{(}\mathbf{x}(J)-\mathbf{x}(0)\big{)}]{\mathrm{d}}\mathbf{x}_{0}\\ &=\lim_{\Delta t\to 0+}\int_{\mathbb{R}^{dn}}\ldots\int_{\mathbb{R}^{dn}}\frac{e^{-\sum_{j=0}^{J-1}\big{(}V(\mathbf{x}(j),\mathbf{X})\Delta t+\frac{|\mathbf{x}({j+1})-\mathbf{x}(j)|^{2}}{2\Delta t}\big{)}}}{(2\pi\Delta t)^{3nJ/2}}\delta\big{(}\mathbf{x}(J)-\mathbf{x}(0)\big{)}\prod_{j-0}^{J}\mathrm{d}\mathbf{x}(j)\\ &=\lim_{\Delta t\to 0+}\int_{\mathbb{R}^{dnJ}}\frac{e^{-\sum_{j=0}^{J-1}\big{(}V(\mathbf{x}(j),\mathbf{X})\Delta t+\frac{|\mathbf{x}({j+1})-\mathbf{x}(j)|^{2}}{2\Delta t}\big{)}}}{(2\pi\Delta t)^{3nJ/2}}{\mathrm{d}}\mathbf{x}(0){\mathrm{d}}\mathbf{x}(1)\ldots{\mathrm{d}}\mathbf{x}(J-1)\,,\\ \end{split}

becomes the Trotter operator splitting, which is $\mathcal{O}((\Delta t)^{2})$ accurate, see [14], due to the boundary condition $\mathbf{x}(J)=\mathbf{x}(0)$ in the right hand side. This integral in $\mathbb{R}^{dnJ}$ is then approximated by the Metropolis or Langevin method applied to the potential

V_{J}(\mathbf{x}):=\sum_{j=0}^{J-1}\Big{(}V\big{(}\mathbf{x}(j),\mathbf{X}\big{)}\Delta t+\frac{|\mathbf{x}({j+1})-\mathbf{x}(j)|^{2}}{2\Delta t}\Big{)}

forming a classical “ring polymer” with “bead” $\mathbf{x}(j)$ , see [15, Chapter 10.2].

2.2. Indistinguishable particles

In the case of $n$ indistinguishable electrons, with the same spin, we have instead the partition function

(2.4)

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{dn}}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\delta(\mathbf{x}_{\beta}-\mathbf{x}_{0}^{\sigma})]{\mathrm{d}}\mathbf{x}_{0}\,,

with the notation $\mathbf{x}^{\sigma}_{t}:=\big{(}\mathbf{x}_{\sigma_{1}}(t),\ldots,\mathbf{x}_{\sigma_{n}}(t)\big{)}$ for $\mathbf{x}_{t}=\big{(}\mathbf{x}_{1}(t),\ldots,\mathbf{x}_{n}(t)\big{)}$ under any permutation $\sigma=(\sigma_{1},\ldots,\sigma_{n})$ of $(1,2,3,\ldots,n)$ and $t\in[0,\beta]$ . The sum $\sum_{\sigma\in\mathbf{S}}$ is over the set of all permutations $\mathbf{S}$ and $\mathbf{x}_{k}(t)\in\mathbb{R}^{d}$ are the coordinates of electron $k$ . This expression of the partition function uses that the potential is invariant with respect to permutation of electron coordinates, i.e $V(\mathbf{x}^{\sigma},\mathbf{X})=V(\mathbf{x},\mathbf{X})$ for all permutations $\sigma$ . The setting with different spin is presented in Appendix C.

The partition function (2.4) can be derived by (2.1) using that the wave functions for the indistinguishable electrons have the antisymmetric representation $\phi(\mathbf{x})=\sum_{\sigma\in\mathbf{S}}{\rm sgn}(\sigma)\tilde{\phi}(\mathbf{x}^{\sigma})/n!$ , for $\tilde{\phi}\in L^{2}(\mathbb{R}^{dn})$ , and consequently

e^{-\beta H_{e}}\phi(\mathbf{x}_{0})=\mathbb{E}[\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\tilde{\phi}(\mathbf{x}_{\beta}^{\sigma})]\,,

which implies (2.4). In the case of bosons the wave functions are symmetric and the trace is replaced by

(2.5)

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{dn}}\sum_{\sigma\in\mathbf{S}}\frac{1}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\delta(\mathbf{x}_{\beta}-\mathbf{x}_{0}^{\sigma})]{\mathrm{d}}\mathbf{x}_{0}\,.

The partition function in (2.4) can also be represented by Brownian bridges as follows. A Brownian bridge process defined by

(2.6)

\mathbf{B}_{t}:=\mathbf{W}_{t}-\frac{t}{\beta}\mathbf{W}_{\beta}\,,\quad 0\leq t\leq\beta\,,

satisfies $\mathbf{B}_{0}=\mathbf{B}_{\beta}=0$ and is independent of $\mathbf{W}_{\beta}$ , since the two Gaussians are uncorrelated

\mathbb{E}[\mathbf{B}_{t}\mathbf{W}_{\beta}^{*}]=\mathbb{E}[(\mathbf{W}_{t}-\frac{t}{\beta}\mathbf{W}_{\beta})\mathbf{W}_{\beta}^{*}]=0\,.

The Wiener path $\mathbf{x}:[0,\beta]\to\mathbb{R}^{3n}$ from $\mathbf{x}_{0}$ to $\mathbf{x}_{\beta}=\mathbf{x}_{0}^{\sigma}$ in (2.4), can be represented by the Brownian bridge

\begin{split}\mathbf{x}_{t}&=\mathbf{x}_{0}+\mathbf{W}_{t}=\mathbf{x}_{0}+\frac{t}{\beta}\mathbf{W}_{\beta}+\mathbf{B}_{t}=(1-\frac{t}{\beta})\mathbf{x}_{0}+\frac{t}{\beta}\mathbf{x}^{\sigma}_{0}+\mathbf{B}_{t}\,.\\ \end{split}

The independence of $\mathbf{B}_{t}=\big{(}\mathbf{B}_{1}(t),\ldots,\mathbf{B}_{n}(t)\big{)}$ and $\mathbf{W}_{\beta}$ implies that the partition function based on the Brownian bridge becomes

(2.7)

\begin{split}\mathrm{Tr}\,(e^{-\beta H_{e}})&=\int_{\mathbb{R}^{3n}}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V({\mathbf{x}}_{t},\mathbf{X}){\mathrm{d}}t}\delta({\mathbf{x}}_{\beta}-\mathbf{x}_{0}^{\sigma})]{\mathrm{d}}\mathbf{x}_{0}\\ &=\int_{\mathbb{R}^{3n}}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\ |\ \mathbf{x}_{\beta}=\mathbf{x}^{\sigma}_{0}]\mathbb{E}[\delta(\mathbf{W}_{\beta}+\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma})]{\mathrm{d}}\mathbf{x}_{0}\\ &=\int_{\mathbb{R}^{3n}}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\ |\ \mathbf{x}_{\beta}=\mathbf{x}^{\sigma}_{0}]\frac{e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}}{(2\pi\beta)^{3n/2}}{\mathrm{d}}\mathbf{x}_{0}\,.\\ \end{split}

The corresponding ring polymer method for fermions following (2.3) becomes

(2.8)

\begin{split}&\mathrm{Tr}\,(e^{-\beta H_{e}})\\ &=\lim_{\Delta t\to 0+}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!(2\pi\Delta t)^{3nJ/2}}\int_{\mathbb{R}^{3n(J+1)}}e^{-V_{J}(\mathbf{x})}\delta\big{(}\mathbf{x}(J)-\mathbf{x}^{\sigma}(0)\big{)}{\mathrm{d}}\mathbf{x}(0){\mathrm{d}}\mathbf{x}(1)\ldots{\mathrm{d}}\mathbf{x}(J)\,.\\ \end{split}

Monte Carlo approximations of the partition functions for bosons and fermions are obtained by sampling both the permutations and the paths, based on the Wiener process. In the case of low temperature, that is $\beta\gg 1$ , the partition function is dominated by the ground state eigenvalue, which typically is substantially lower for bosons than for fermions. Consequently, the value of the trace is lower for fermions than for bosons, while the variance in the Monte Carlo methods for the two remains similar unless the cancellations for fermions can be removed by rewriting the fermions partition function as a sum with few cancelling terms. The fermion sign problem refers to a much larger computational work to approximate the fermion partition function by Monte Carlo samples of permutations and paths, compared to the boson partition function, in the case the sign cancellation is not explicitly included, cf. [8, 11, 5]. The computationally costly sampling of all permutations has been improved by using the worm algorithm [3].

The main new mathematical idea here is to write the fermion partition function (2.7) as an integral of a determinant based on Brownian bridge paths. We thereby reduce the sign problem regarding the complexity with respect to the number of particles. The sign problem related to the increasing computational work with large $\beta$ , however, remains.

Fermion determinants have been formulated before, both for quantum lattice models using hopping between discrete lattice points in a second quantization setting, see the determinant quantum Monte Carlo method [2, 23], and for path integrals in the first quantization setting, based on the ring polymer approximations of (2.8), see [31, 25] and (3.1). Our formulation also uses the path integral in the first quantization formulation (2.7) but not the ring polymer approximation, instead we use a different determinant representation based on explicit generation of Brownian bridge paths from the Wiener process. The first quantization formulation here, with Brownian bridge particle paths in $\mathbb{R}^{3n}$ , has the two advantages: (1) simple statistical independent sampling of Wiener processes as compared to Metropolis or Langevin sampling of ring polymers in high dimension, and (2) the computational work, for $n$ particles requiring $n$ classical paths, is proportional to $\mathcal{O}(n^{3})$ . The drawback is that in our formulation two body interactions cannot be represented exactly. Nonetheless, we provide a computable error indicator based on perturbation analysis.

For fermion particles in dimension one, $d=1$ , the work [9, 13] shows that reflected Wiener processes can be used to avoid the fermion sign problem in (2.4). Thereby the computational complexity becomes roughly $\mathcal{O}(\epsilon^{-3})$ , as for distinguishable particles, as a function of the expected approximation error $\epsilon$ . However, for fermion particles in higher dimension, $d>1$ , it seems that the literature only provides computational methods to approximate (2.4) with a computational complexity that grows exponentially with the number of particles $n$ , cf. [11], although there are different approximation alternatives, e.g., [30, 9, 25, 19], that improve the computational work compared to a straightforward Monte Carlo approximation of (2.4).

Simulations of bosons in (2.5) do not experience the sign problem, nonetheless, they are also hard. For instance, the related problem to determine a matrix permanent is NP-hard, see [18].

The fixed node method [1, 10, 7] avoids the fermion sign problem for electron ground states, related to the formulation in $d=1$ , but requires the nodal set of the wave function, which itself is challenging to determine and consequently causes computational bias that is difficult to estimate.

Our study relies on the electron operator, $H_{e}=-\frac{1}{2}\Delta_{\mathbf{x}}+V$ , and uses specific properties of both the Laplacian and the Coulomb interaction potentials for electrons and nuclei.

3. The fermion partition functions based on determinants

Takahashi and Imada [31] applied the fermion antisymmetrization to the density matrix approximation for each time step, namely to $e^{-(\frac{|\mathbf{x}(j+1)-\mathbf{x}(j)|^{2}}{2\Delta t}+V(\mathbf{x}(j),\mathbf{X})\Delta t)}$ in (2.3), to obtain a formulation of the partition function based on products of determinants as follows

(3.1)

\begin{split}&\mathrm{Tr}\,(e^{-\beta H_{e}})\\ &=\lim_{\Delta t\to 0+}\int_{\mathbb{R}^{3nJ}}\prod_{j=0}^{J-1}\Big{(}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!(2\pi\Delta t)^{3n/2}}e^{-(\frac{|\mathbf{x}^{\sigma}(j+1)-\mathbf{x}(j)|^{2}}{2\Delta t}+V(\mathbf{x}(j),\mathbf{X})\Delta t)}\Big{)}{\mathrm{d}}\mathbf{x}(0)\ldots{\mathrm{d}}\mathbf{x}(J-1)\\ &=\lim_{\Delta t\to 0+}\int_{\mathbb{R}^{3nJ}}\prod_{j=0}^{J-1}\Big{(}\frac{e^{-V(\mathbf{x}(j),\mathbf{X})\Delta t}}{n!(2\pi\Delta t)^{3n/2}}{\rm det}\mathcal{R}\big{(}\mathbf{x}(j),\mathbf{x}(j+1)\big{)}\Big{)}{\mathrm{d}}\mathbf{x}(0)\ldots{\mathrm{d}}\mathbf{x}(J-1)\,,\\ \end{split}

where $\mathbf{x}(J)=\mathbf{x}(0)$ and the $n\times n$ matrix $\mathcal{R}$ has the components

\mathcal{R}_{\ell k}\big{(}\mathbf{x}(j),\mathbf{x}(j+1)\big{)}:=e^{|\mathbf{x}_{k}(j+1)-\mathbf{x}_{\ell}(j)|^{2}/(2\Delta t)}\,.

A determinant can be computed with $\mathcal{O}(n^{3})$ work. Therefore the formulation avoids the $\mathcal{O}(n!)$ computational complexity to sum over all permutations. However since the determinants have varying sign, for particles in dimension two and higher, the required Monte Carlo sampler suffers from the fermion sign problem, which leads to a large statistical variance [31, 25].

In this section we present a different determinant formulation which is based on (2.7) where the Monte Carlo method directly samples independent Brownian bridge paths, with the aim to reduce the fermion sign problem.

The formulation of the partition function

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{3n}}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\ |\ \mathbf{x}_{\beta}=\mathbf{x}^{\sigma}_{0}]\frac{e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}}{(2\pi\beta)^{3n/2}}{\mathrm{d}}\mathbf{x}_{0}

using the Brownian bridge $\mathbf{x}_{t}=\big{(}1-\frac{t}{\beta}\big{)}\mathbf{x}_{0}+\frac{t}{\beta}\mathbf{x}_{0}^{\sigma}+\mathbf{B}_{t}$ can be represented with a determinant as follows. Assume first that the potential is particle wise separable, so that

(3.2)

V(\mathbf{x},\mathbf{X})=\sum_{k=1}^{n}\tilde{V}_{k}(\mathbf{x}_{k},\mathbf{X}),

then

H_{e}=-\frac{1}{2}\Delta_{\mathbf{x}}+V=\sum_{k=1}^{n}\big{(}-\frac{1}{2}\Delta_{\mathbf{x}_{k}}+\tilde{V}_{k}(\mathbf{x}_{k},\mathbf{X})\big{)}

is separable. Define the $n\times n$ matrix

\mathcal{W}_{k\ell}\big{(}\mathbf{x}\big{)}:=e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}e^{-\int_{0}^{\beta}\tilde{V}_{k}(\mathbf{B}_{k}(t)+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\ell}(0),\mathbf{X}){\mathrm{d}}t}\,.

Its determinant yields the path integral

\begin{split}&{\mathrm{det}}\big{(}\mathcal{W}(\mathbf{x})\big{)}=\sum_{\sigma\in\mathcal{S}}{\rm sgn}(\sigma)\mathcal{W}_{1\sigma_{1}}\mathcal{W}_{2\sigma_{2}}\ldots\mathcal{W}_{n\sigma_{n}}\\ &=\sum_{\sigma\in\mathcal{S}}{\rm sgn}(\sigma)e^{-\sum_{k=1}^{n}|\mathbf{x}_{k}(0)-\mathbf{x}_{\sigma_{k}}(0)|^{2}/(2\beta)}e^{-\int_{0}^{\beta}\sum_{k=1}^{n}\tilde{V}_{k}(\mathbf{B}_{k}(t)+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\sigma_{k}}(0),\mathbf{X}){\mathrm{d}}t}\\ &=\sum_{\sigma\in\mathcal{S}}{\rm sgn}(\sigma)e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}e^{-\int_{0}^{\beta}V(\mathbf{B}(t)+(1-\frac{t}{\beta})\mathbf{x}(0)+\frac{t}{\beta}\mathbf{x}^{\sigma}(0),\mathbf{X}){\mathrm{d}}t}\,,\\ \end{split}

and we obtain the representation

(3.3)

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{3n}}\mathbb{E}[\frac{{\rm det}\big{(}\mathcal{W}(\mathbf{x})\big{)}}{n!(2\pi\beta)^{3n/2}}]{\mathrm{d}}\mathbf{x}_{0}\,.

The path integral can be approximated by the trapezoidal method as follows. Make a partition $t_{m}=m\Delta t$ , for $m=0,\ldots M$ , with $\Delta t=\beta/M$ and let

\begin{split}\mathcal{\overline{W}}_{k\ell}\big{(}\mathbf{x}\big{)}&:=e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}e^{-\sum_{m=1}^{M-1}\tilde{V}_{k}(\mathbf{B}_{k}(t_{m})+(1-\frac{t_{m}}{\beta})\mathbf{x}_{k}(0)+\frac{t_{m}}{\beta}\mathbf{x}_{\ell}(0))\Delta t}\times\\ &\quad\times e^{-\tilde{V}(\mathbf{B}_{k}(0)+\mathbf{x}_{k}(0))\Delta t/2-\tilde{V}_{k}(\mathbf{B}_{k}(0)+\mathbf{x}_{\ell}(0))\Delta t/2}\,.\end{split}

Then we have the computable second order accurate approximation

(3.4)

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{3n}}\mathbb{E}[\frac{{\rm det}\big{(}\mathcal{\overline{W}}(\mathbf{x})\big{)}}{n!(2\pi\beta)^{3n/2}}]{\mathrm{d}}\mathbf{x}_{0}+\mathcal{O}((\Delta t)^{2})\,.

Such second order accuracy of path integral approximations is proved in [14] assuming that the functional, here $\mathbf{x}\mapsto e^{-\int_{0}^{\beta}V(\mathbf{x}_{t}){\mathrm{d}}t}$ , is six times Fréchet differentiable. The determinant for an $n\times n$ matrix can be determined with $\mathcal{O}(n^{3})$ number of operations using the $LU$ -factorization in contrast to the $\mathcal{O}(n!)$ computational work to sum over all permutations.

The Hamiltonian for molecular systems has a potential based on Coulomb interactions of the electrons and nuclei, see [6, 22, 24],

(3.5)

\begin{split}H_{e}(\mathbf{x},\mathbf{X})&=-\frac{1}{2}\Delta_{\mathbf{x}}+V(\mathbf{x},\mathbf{X})\,,\\ V(\mathbf{x},\mathbf{X})&=\sum_{k=1}^{n}\sum_{k<k^{\prime}}\frac{1}{|\mathbf{x}_{k}-\mathbf{x}_{k^{\prime}}|}+\sum_{j=1}^{N}\sum_{j<i}\frac{Z_{j}Z_{i}}{|\mathbf{X}_{j}-\mathbf{X}_{i}|}-\sum_{k=1}^{n}\sum_{j=1}^{N}\frac{Z_{j}}{|\mathbf{X}_{j}-\mathbf{x}_{k}|}\,,\end{split}

where $Z_{i}$ denotes the charge of the $i$ th nucleus. The coordinates $\mathbf{x}_{k}\in\mathbb{R}^{3}$ and $\mathbf{X}_{j}\in\mathbb{R}^{3}$ are for electron $k$ and nucleus $j$ . In this case the terms depending on $\mathbf{x}$ can still be written as a sum over fermion particles

\sum_{k=1}^{n}(\sum_{k^{\prime}\neq k}\frac{1}{2|\mathbf{x}_{k}-\mathbf{x}_{k^{\prime}}|}-\sum_{j=1}^{N}\frac{Z_{j}}{|\mathbf{x}_{k}-\mathbf{X}_{j}|})=:\sum_{k=1}^{n}V_{k}(\mathbf{x},\mathbf{X})

to form the potential energy per particle, $V_{k}$ , but the electron repulsion sum is not separable particle wise and the repulsion, $\sum_{k^{\prime}\neq k}\frac{1}{2|\mathbf{x}_{k}-\mathbf{x}_{k^{\prime}}|}$ , is instead a sum including all row coordinates $\mathbf{x}_{k^{\prime}}$ .

Let

\mathbf{x}_{k}^{\sigma_{k}}(t):=\mathbf{B}_{k}(t)+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\sigma_{k}}(0)\,,

then we have

(3.6)

\begin{split}\mathrm{Tr}\,(e^{-\beta H_{e}})&=\int_{\mathbb{R}^{3n}}\mathbb{E}\big{[}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!(2\pi\beta)^{3n/2}}e^{-\sum_{k=1}^{n}|\mathbf{x}_{k}(0)-\mathbf{x}_{\sigma_{k}}(0)|^{2}/(2\beta)}\times\\ &\quad\times e^{-\frac{\beta}{2}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}Z_{i}Z_{j}/|\mathbf{X}_{i}-\mathbf{X}_{j}|}\times\\ &\quad\times e^{\sum_{k=1}^{n}\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}\times\\ &\quad\times e^{-\frac{1}{2}\sum_{k=1}^{n}\int_{0}^{\beta}\sum_{k^{\prime}=1,k^{\prime}\neq k}^{n}|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{k^{\prime}}^{\sigma_{k^{\prime}}}(t)|^{-1}{\mathrm{d}}t}\big{]}{\mathrm{d}}\mathbf{x}(0)\,,\end{split}

which cannot be written as a determinant of an $n\times n$ matrix, for $n>2$ , due to the electron repulsion. Instead it is a sum over permutations for a tensor. Sums of permutations of tensors are often computationally demanding, e.g. $NP$ -hard see [18]. An exception is the computation of the determinant for an $n\times n$ matrix which can be determined with $\mathcal{O}(n^{3})$ operations. The next step is therefore to formulate approximations to the right hand side in (3.6), avoiding to sum over all permutations.

In the special case with two particles, $n=2$ , the electron repulsion in (3.6) can in fact be written as a determinant, since

(3.7)

\sigma_{k^{\prime}}=\left\{\begin{array}[]{cc}2&\mbox{ if }\sigma_{k}=1\,,\\ 1&\mbox{ if }\sigma_{k}=2\,.\\ \end{array}\right.

A natural approximation of (3.6) is to replace $1/|\mathbf{x}_{k}^{\sigma_{k}}-\mathbf{x}_{k^{\prime}}^{\sigma_{k^{\prime}}}|$ by $1/|\mathbf{x}_{k}^{\sigma_{k}}-\mathbf{x}_{k^{\prime}}^{\nu_{k^{\prime}}}|$ where

\nu_{k^{\prime}}=\left\{\begin{array}[]{ll}k^{\prime},&\mbox{ if }k^{\prime}\neq\sigma_{k}\,,\\ k,&\mbox{ if }k^{\prime}=\sigma_{k}\,,\end{array}\right.

which for $n=2$ becomes the exact form (3.7) and for $n>2$ treats the remaining particles with coordinates $\mathbf{x}_{k^{\prime}}$ as distinguishable as follows. Define the approximation

(3.8)

\begin{split}\mathcal{Z}_{\nu}&:=\int_{\mathbb{R}^{3n}}\mathbb{E}\big{[}\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!(2\pi\beta)^{3n/2}}e^{-\sum_{k=1}^{n}|\mathbf{x}_{k}(0)-\mathbf{x}_{\sigma_{k}}(0)|^{2}/(2\beta)}\times\\ &\quad\times e^{-\frac{\beta}{2}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}Z_{i}Z_{j}/|\mathbf{X}_{i}-\mathbf{X}_{j}|}\times\\ &\quad\times e^{\sum_{k=1}^{n}\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}\times\\ &\quad\times e^{-\frac{1}{2}\sum_{k=1}^{n}\int_{0}^{\beta}\sum_{j=1,j\neq k}^{n}|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|^{-1}{\mathrm{d}}t}\big{]}{\mathrm{d}}\mathbf{x}(0)\,,\end{split}

which can be formulated by a determinant: let

\begin{split}\mathcal{W}_{k\ell}(\mathbf{x})&:=e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}e^{\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\ell}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j=1,j\neq k}^{n}|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|^{-1}{\mathrm{d}}t}\,,\end{split}

where

\nu_{j}=\left\{\begin{array}[]{ll}j,&\mbox{ if }j\neq\ell\,,\\ k,&\mbox{ if }j=\ell\,,\end{array}\right.\,

then we have

(3.9)

\mathcal{Z}_{\nu}=e^{-\frac{\beta}{2}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}Z_{i}Z_{j}/|\mathbf{X}_{i}-\mathbf{X}_{j}|}\int_{\mathbb{R}^{3n}}\mathbb{E}\big{[}\frac{{\rm det}\big{(}\mathcal{W}(\mathbf{x})\big{)}}{n!(2\pi\beta)^{3n/2}}\big{]}{\mathrm{d}}\mathbf{x}(0)\,,

and $\mathcal{Z}_{\nu}=\mathrm{Tr}\,(e^{-\beta H_{e}})$ for $n=2$ . We also note that for $\sigma$ equal to the identity permutation or a permutation with only one transposition it holds that

|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\sigma_{j}}(t)|=|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|\,.

For small values of $\beta$ such permutations will dominate in (3.8) and (3.6), due to the factors $e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\sigma_{k}}(0)|^{2}/(2\beta)}$ , and consequently (3.8) is a consistent approximation of (3.6) as $\beta\to 0+$ .

3.1. Motivation for the determinant approximation (3.8)

The aim of this section is to motivate that for a given particle $k$ the interaction

e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j\neq k}|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\sigma_{j}}(t)|^{-1}{\mathrm{d}}t}

in (3.6), where $\mathbf{x}_{i}^{\sigma_{i}}(t)=\mathbf{B}_{i}(t)+(1-\frac{t}{\beta})\mathbf{x}_{i}(0)+\frac{t}{\beta}\mathbf{x}_{\sigma_{i}}(0)$ , can be related to the interaction in the determinant formulation

e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j\neq k}|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|^{-1}{\mathrm{d}}t}\,.

In particular we study the marginal distribution of $\mathbf{x}_{\sigma_{j}}(0)$ for Brownian bridges from $\mathbf{x}(0)$ to $\mathbf{x}^{\sigma}(0)$ . The mean of this marginal distribution motivates the choice $\nu_{j}$ . This section shows that the variances of the marginals lead to the perturbed interaction

(3.10)

\mathbb{E}_{\xi}[e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j\neq k}|\mathbf{x}_{k}^{\sigma_{k}}(t)-(\mathbf{x}_{j}^{\nu_{j}}(t)+\frac{t}{\beta}\sqrt{\frac{\beta}{c_{*}}}\,\xi_{j})|^{-1}{\mathrm{d}}t}]\,,

which we use for sensitivity analysis of the determinant formulation. Here the expected value is with respect to only $\xi_{j}$ , for $j\neq k$ , which are independent standard normal distributed random variables in $\mathbb{R}^{3}$ and $c_{*}\approx 2$ is a parameter related to the logarithm of the permutation cycle length. In particular, the numerical sensitivity experiments in Table 4.6 show that the perturbed mean-field approximation based on (3.10) has a relative error on the same order of accuracy as $\bar{h}_{\nu}$ and thereby provides a computational error indicator to roughly estimate the accuracy of $\bar{h}_{\nu}$ , without using a reference value.

3.1.1. Marginals of $\mathbf{x}_{\sigma_{j}}(0)$

If $\sigma={\rm I}$ we have $\mathbf{x}_{j}^{\sigma_{j}}=\mathbf{x}_{j}^{\nu_{j}}$ and for small $\beta$ the permutation $\sigma={\rm I}$ dominates due to the factor $e^{-|\mathbf{x}(0)-\mathbf{x}^{\sigma}(0)|^{2}/(2\beta)}$ in (3.6). The case $\sigma={\rm I}$ implies to have only one-cycle permutations. In the case of a two cycle for particle $j$ , i.e. $\sigma_{\sigma_{j}}=j$ , the weight involving $\mathbf{x}_{\sigma_{j}}(0)$ becomes

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{j}}(0)-\mathbf{x}_{\sigma_{\sigma_{j}}}(0)|^{2}/(2\beta)}=e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/\beta}\,.

For a three-cycle, $\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)=\mathbf{x}_{j}(0)$ , we have

\begin{split}&e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{j}}(0)-\mathbf{x}_{\sigma_{\sigma_{j}}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{\sigma_{j}}}(0)-\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)|^{2}/(2\beta)}\\ &=e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{\sigma_{j}}}(0)-\frac{\mathbf{x}_{j}(0)+\mathbf{x}_{\sigma_{j}}(0)}{2\beta}|^{2}/\beta}e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(4\beta)}\,,\end{split}

using that for any $y,y_{a}$ and $y_{b}$ in $\mathbb{R}^{3}$

e^{-|y-y_{a}|^{2}/\gamma_{a}}e^{-|y-y_{b}|^{2}/\gamma_{b}}=e^{-|y-\frac{\gamma_{b}y_{a}+\gamma_{a}y_{b}}{\gamma_{a}\gamma_{b}}|^{2}\frac{\gamma_{a}+\gamma_{b}}{\gamma_{a}\gamma_{b}}}e^{-|y_{a}-y_{b}|^{2}/(\gamma_{a}+\gamma_{b})}\,.

Integration with respect to $\mathbf{x}_{\sigma_{\sigma_{j}}(0)}$ yields the non normalized marginal distribution

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}(\frac{1}{2\beta}+\frac{1}{4\beta})}=e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\frac{3}{4\beta}}\,.

For a four-cycle with $\mathbf{x}_{\sigma_{\sigma_{\sigma_{\sigma_{j}}}}}(0)=\mathbf{x}_{j}(0)$ we obtain similarly the weight

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{j}}(0)-\mathbf{x}_{\sigma_{\sigma_{j}}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{\sigma_{j}}}(0)-\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)-\mathbf{x}_{\sigma_{\sigma_{\sigma_{\sigma_{j}}}}}(0)|^{2}/(2\beta)}

which by integration with respect to $\mathbf{x}_{\sigma_{\sigma_{j}}}(0)$ and $\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)$ determines the marginal distribution proportional to

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}(\frac{1}{2\beta}+\frac{1}{4\beta}+\frac{1}{4\beta+2\beta})}=e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\frac{11}{12\beta}}\,.

For an $m$ -cycle, $m\geq 3$ , we obtain analogously the marginal distribution proportional to

(3.11)

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\sum_{m^{\prime}=1}^{m-1}\frac{1}{2m^{\prime}\beta}}\,,

where

(3.12)

\sum_{m^{\prime}=1}^{m-1}\frac{1}{2m^{\prime}\beta}=:\frac{c_{m}}{2\beta}\approx\frac{c_{*}}{2\beta}\,,

and $c_{2}=2,c_{3}=3/2,c_{4}=11/6,c_{5}=25/12\,.$ The approximation $c_{*}=2$ for small $\beta$ is further motivated in Section 3.1.2.

In the determinant formulation (3.9)

\mathcal{W}_{k\ell}=e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}e^{\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\ell}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j=1,j\neq k}^{n}|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|^{-1}{\mathrm{d}}t}\,,

we do not have access to the permutations. Therefore we approximate $\mathbf{x}_{\sigma_{j}}(0)$ in the interaction $|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\sigma_{j}}(t)|^{-1}$ by the mean $\mathbf{x}_{j}(0)$ of the normal distributions $e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\sum_{m^{\prime}=1}^{m-1}\frac{1}{2m^{\prime}\beta}}$ , with respect to the variable $\mathbf{x}_{\sigma_{j}}(0)$ , related to $m$ -cycles. If $j=\ell$ and $\ell\neq k$ we cannot use $\nu_{j}=\ell$ since then $|\mathbf{x}_{k}^{\ell}(\beta)-\mathbf{x}_{\ell}^{\ell}(\beta)|=0$ and consequently we would obtain $\mathcal{W}_{k\ell}=0$ for all $\ell$ in the numerical approximations. Instead we let $\nu_{\ell}=k$ in order for $\nu_{j},\,j\neq k$ , to all be different.

3.1.2. The perturbation formulation (3.10)

To motivate the formulation (3.10), we introduce the notation

\begin{split}\mathrm{Tr}\,(e^{-\beta H_{e}})&=\sum_{\sigma\in\mathbf{S}}\frac{{\rm sgn}(\sigma)}{n!(2\pi\beta)^{3n/2}}\int_{\mathbb{R}^{3n}}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\ |\ \mathbf{x}_{\beta}=\mathbf{x}^{\sigma}_{0}]e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{0}\\ &=:\sum_{\sigma\in\mathbf{S}}{\rm sgn}(\sigma)\int_{\mathbb{R}^{3n}}f_{\sigma}(\mathbf{x}_{0})e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{0}\,.\\ \end{split}

Assume that for a permutation $\sigma$ the coordinate $\mathbf{x}_{j}(0)$ is in an $m$ -cycle consisting of the coordinates

\big{(}\mathbf{x}_{j}(0),\mathbf{x}_{\sigma_{j}}(0),\ldots,\mathbf{x}_{\sigma_{\ldots\sigma_{j}}}(0)\big{)}=:y_{\sigma}\in\mathbb{R}^{3m}

and let the particle coordinates that are not in this $m$ -cycle be denoted by $y^{\perp}_{\sigma}\in\mathbb{R}^{3(n-m)}$ . Define

f_{\sigma}\big{(}\mathbf{x}_{1}(0),\mathbf{x}_{2}(0),\ldots,\mathbf{x}_{n}(0)\big{)}=:f_{\sigma,j}(y_{\sigma},y^{\perp}_{\sigma})\,.

We have by (3.11)

\begin{split}&\int_{\mathbb{R}^{3n}}f_{\sigma}(\mathbf{x}_{0})e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{0}\\ &=\int_{\mathbb{R}^{3n}}f_{\sigma,j}(y_{\sigma},y^{\perp}_{\sigma})e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{j}}(0)-\mathbf{x}_{\sigma_{\sigma_{j}}}(0)|^{2}/(2\beta)}e^{-|\mathbf{x}_{\sigma_{\sigma_{j}}}(0)-\mathbf{x}_{\sigma_{\sigma_{\sigma_{j}}}}(0)|^{2}/(2\beta)}\ldots\times\\ &\qquad\times e^{-|y^{\perp}-(y^{\perp})^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{j}(0){\mathrm{d}}\mathbf{x}_{\sigma_{j}}(0)\ldots{\mathrm{d}}\mathbf{x}_{\sigma_{\ldots\sigma_{j}}}(0){\mathrm{d}}y^{\perp}_{\sigma}\\ &=C_{mn}\int_{\mathbb{R}^{3(n+1-m)}}f_{\sigma,j}\Big{(}\mathbf{x}_{j}(0),\mathbf{x}_{\sigma_{j}}(0),\bar{\mathbf{x}}_{\sigma_{\sigma_{j}}}(0),\bar{\mathbf{x}}_{\sigma_{\sigma_{\sigma_{j}}}}(0),\ldots,\bar{\mathbf{x}}_{\sigma_{\ldots\sigma_{j}}}(0),y^{\perp}_{\sigma}\Big{)}\times\\ &\qquad\times e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\sum_{m^{\prime}=1}^{m-1}{\frac{1}{2m^{\prime}\beta}}}e^{-|y^{\perp}-(y^{\perp})^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{j}(0){\mathrm{d}}\mathbf{x}_{\sigma_{j}}(0){\mathrm{d}}y^{\perp}_{\sigma}\,,\end{split}

where $\big{(}\bar{\mathbf{x}}_{\sigma_{\sigma_{j}}}(0),\bar{\mathbf{x}}_{\sigma_{\sigma_{\sigma_{j}}}}(0),\ldots,\bar{\mathbf{x}}_{\sigma_{\ldots\sigma_{j}}}(0)\big{)}$ and the constant $C_{mn}$ are defined by the integration with respect to $\big{(}\mathbf{x}_{\sigma_{\sigma_{j}}}(0),\ldots,\mathbf{x}_{\sigma_{\ldots\sigma_{j}}}(0)\big{)}\in\mathbb{R}^{3(m-2)}\,,$ using the mean value theorem. The mean value theorem also yields

(3.13)

\begin{split}&\int_{\mathbb{R}^{3n}}f_{\sigma}(\mathbf{x}_{0})\frac{e^{-|\mathbf{x}_{0}-\mathbf{x}_{0}^{\sigma}|^{2}/(2\beta)}}{(2\pi\beta)^{3n/2}}{\mathrm{d}}\mathbf{x}_{0}\\ &=C^{\prime}_{mn}\int_{\mathbb{R}^{3(n+1-m)}}f_{\sigma,j}\Big{(}\mathbf{x}_{j}(0),\bar{\mathbf{x}}_{\sigma_{j}}(0),\bar{\mathbf{x}}^{\prime}_{\sigma_{\sigma_{j}}}(0),\bar{\mathbf{x}}^{\prime}_{\sigma_{\sigma_{\sigma_{j}}}}(0),\ldots,\bar{\mathbf{x}}^{\prime}_{\sigma_{\ldots\sigma_{j}}}(0),y^{\perp}_{\sigma}\Big{)}\times\\ &\qquad\times e^{-|y^{\perp}-(y^{\perp})^{\sigma}|^{2}/(2\beta)}{\mathrm{d}}\mathbf{x}_{j}(0){\mathrm{d}}y^{\perp}_{\sigma}\,,\end{split}

by integration with respect to the Gaussian measure $e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\sum_{m^{\prime}=1}^{m-1}{\frac{1}{2m^{\prime}\beta}}}{\mathrm{d}}\mathbf{x}_{\sigma_{j}}(0)$ , for another constant $C^{\prime}_{mn}$ and new values $(\bar{\mathbf{x}}^{\prime}_{\sigma_{\sigma_{j}}}(0),\bar{\mathbf{x}}^{\prime}_{\sigma_{\sigma_{\sigma_{j}}}}(0)),\ldots,\bar{\mathbf{x}}^{\prime}_{\sigma_{\ldots\sigma_{j}}}(0))$ obtained from the mean value theorem and the dependence on $\mathbf{x}_{\sigma_{j}}(0)$ .

The work [12] shows that, provided the temperature is high enough, $\beta\lesssim 1$ , the probability of cycle length $m$ with high probability decays exponentially for the ideal Fermi gas, i.e. for non interacting electrons, and for the uniform electron gas, i.e. for interacting electrons. On the other hand, in the case of low temperature, $\beta=8$ , Figure 9 in [12] shows that the probability is almost constant with respect to the cycle length. In our study focused on high temperature, we therefore approximate the sum in (3.11) and (3.12) as

e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\sum_{m^{\prime}=1}^{m-1}\frac{1}{2m^{\prime}\beta}}\approx e^{-|\mathbf{x}_{j}(0)-\mathbf{x}_{\sigma_{j}}(0)|^{2}\frac{c_{*}}{2\beta}}\,,

where $c_{*}=2$ . The dependence on $\sigma$ in (3.11) and (3.12), through the cycle length $m$ , is thereby avoided so that we can apply a determinant formulation. The normal distribution of $\mathbf{x}_{\sigma_{j}}(0)$ obtained in (3.11) and (3.13) provides a method to study perturbations of the determinant formulation (3.9), where the size of the perturbations are related to the unavailable interactions $|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{\sigma_{j}}(t)|$ for the determinant formulation. Therefore we use in (3.10) the $\sigma$ independent approximation

(3.14)

\bar{\mathbf{x}}_{\sigma_{j}}(0)\approx\mathbf{x}_{j}(0)+\sqrt{\frac{\beta}{c_{*}}}\,\xi_{j}

where $\xi_{j}\in\mathbb{R}^{3}$ are independent standard normal random variables. A corresponding rigorous analysis of the perturbations would require to more precisely determine the approximation error in (3.14) using the independent samples $\xi_{j}$ . The perturbation approximation (3.14) motivates to use the following error indicator for the determinant formulation

(3.15)

\begin{split}&|\mathrm{Tr}\,(e^{-\beta H_{e}})-\mathcal{Z}_{\nu}|\\ &=C\Big{|}\sum_{\sigma\in\mathcal{S}}{\rm sgn}(\sigma)\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\prod_{k=1}^{n}\Big{(}e^{\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}\times\\ &\qquad\times e^{-\sum_{j\neq k}\frac{1}{2}\int_{0}^{\beta}|\mathbf{x}_{k}^{\sigma_{k}}(t)-\mathbf{x}_{j}^{j}(t)-(\mathbf{x}_{j}^{\sigma_{j}}(t)-\mathbf{x}_{j}^{j}(t))|^{-1}{{\mathrm{d}}t}}e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\sigma_{k}}(0)|^{2}/(2\beta)}\Big{)}\Big{]}{\mathrm{d}}\mathbf{x}_{0}\\ &\quad-\int_{\mathbb{R}^{3n}}\mathbb{E}[{\rm det}\big{(}\mathcal{W}(\mathbf{x})\big{)}]{\mathrm{d}}\mathbf{x}_{0}\Big{|}\\ &\approx C|\int_{\mathbb{R}^{3n}}\mathbb{E}[{\rm det}\big{(}\widetilde{\mathcal{W}}(\mathbf{x})\big{)}]{\mathrm{d}}\mathbf{x}_{0}-\int_{\mathbb{R}^{3n}}\mathbb{E}[{\rm det}\big{(}\mathcal{W}(\mathbf{x})\big{)}]{\mathrm{d}}\mathbf{x}_{0}|\end{split}

where

(3.16)

\begin{split}\widetilde{\mathcal{W}}_{k\ell}(\mathbf{x}_{0})&=e^{\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\ell}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}\times\\ &\qquad\times\mathbb{E}_{\xi}[e^{-\sum_{j\neq k}\frac{1}{2}\int_{0}^{\beta}{|\mathbf{x}_{k}^{\ell}(t)-(\mathbf{x}_{j}^{\nu_{j}}(t)+\frac{t}{\beta}\sqrt{\frac{\beta}{c_{*}}}\xi_{j})|^{-1}}{{\mathrm{d}}t}}]e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}\,,\\ C&:=\frac{e^{-\frac{\beta}{2}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}Z_{i}Z_{j}/|\mathbf{X}_{i}-\mathbf{X}_{j}|}}{n!(2\pi\beta)^{3n/2}}\,.\end{split}

Here the expected value of ${\rm det}(\widetilde{\mathcal{W}})$ is with respect to $\mathbf{B}$ and the expeded value $\mathbb{E}_{\xi}$ is with respect to $\xi$ only.

Remark 3.1 (Bias in mean-field approximations).

In the case of separable potentials the determinant mean-field formulation (3.9) and the Hartree-Fock mean-field method for electron ground state approximations both have no bias. The determinant formulation also has no bias in the case of interacting distinguishable particles, which is not the case for the Hartree method. Consequently the determinant formulation is asymptotically accurate in the two important settings of separable fermion potentials and for interacting distinguishable particles.

4. Numerical experiments

In this section, we apply the path integral with Feynman-Kac formulation (3.4) and (3.9), to approximate the quantum statistical partition function and the mean-field energy for systems consisting of a few indistinguishable electrons. These systems operate at various inverse temperatures $\beta$ and involve two different potential functions denoted as $V(\mathbf{x},\,\mathbf{X})$ , namely:


(V1)	$\displaystyle V(\mathbf{x},\mathbf{X})$	$\displaystyle=\frac{1}{2}\|\mathbf{x}\|^{2},\textup{ Section }\ref{Subsection:Fermi_gas},$
(V2)	$\displaystyle V(\mathbf{x},\mathbf{X})$	$\displaystyle=\frac{1}{2}\|\mathbf{x}\|^{2}+\frac{1}{2}\sum_{k}\sum_{\ell\neq k}\frac{\lambda}{\|x_{k}-x_{\ell}\|},\textup{ Section }\ref{Subsection_Coulomb_repulsion_plus_HO}.$

The derived partition function approximations facilitate the computation of corresponding mean-field energy values. To validate our approach, we compare our results for the two test cases, V1 and V2, with three different benchmarks: exact solutions obtained for particle-wise separable potentials, reference data in [11], and the formula (2.7) employing all the permutations. Furthermore, we implement the perturbation formulation in Section 3.1.2 to calculate an error indicator for our numerical results. The MATLAB code for implementing the algorithm based on (3.4) and (3.9) is openly available through our GitHub repository [34].

4.1. Non-interacting Fermi gas in a confining harmonic trap

In $d$ -dimensional space with a governing separable external potential, we study a system comprising $n$ indistinguishable non-interacting fermions, known as a Fermi gas. Our investigation involves a test case where the Fermi gas is confined by an external harmonic trap defined by $V_{\mathrm{HO}}(\mathbf{x})=\frac{1}{2}|\mathbf{x}|^{2}$ . To determine the partition function’s exact value at inverse temperature $\beta=1\,/\,(k_{B}T)$ , we use a recursive formula, given by

(4.2)

Z_{n}(\beta)=\frac{1}{n}\,\sum_{k=1}^{n}\,(-1)^{k-1}\,Z_{1}(k\beta)\,Z_{n-k}(\beta)\,,\mbox{ with }Z_{0}(\beta)=1\,,

where

(4.3)

Z_{1}(\beta)=\sum_{j=1}^{\infty}\,e^{-\beta\,\epsilon_{j}}\,,

represents the single-particle partition function. Here, $\epsilon_{j}$ denotes the energy levels of each eigenstate corresponding to $V_{\mathrm{HO}}(\mathbf{x})$ , see [4, 28]. We compute the exact partition function using (4.2) as the reference value, which serves to test the convergence of the proposed path integral Monte Carlo method (3.4).

Considering the simple case of a one-dimensional external harmonic oscillator potential, the single-particle energy eigenvalues in Hartree atomic units are expressed as $\eta_{j}=\frac{1}{2}+j$ , where $j=0,1,2,\dots$ . Using these eigenvalues, we can directly deduce the single-particle partition function $Z_{1}(\beta)$ for a three-dimensional harmonic oscillator trap:

(4.4)

Z_{1}(\beta)=\sum_{j,k,\ell=0}^{\infty}e^{-\beta(\eta_{j}+\eta_{k}+\eta_{\ell})}=\Big{(}\sum_{j=0}^{\infty}e^{-\beta\eta_{j}}\Big{)}^{3}=\Big{(}\frac{\exp{(-\frac{\beta}{2})}}{1-\exp{(-\beta)}}\Big{)}^{3}\,.

By combining the recursive formula (4.2) with (4.4), we can compute the exact partition function $Z_{n}(\beta)$ . In contrast, using (3.9), we can approximate the partition function for the test case without nuclei by

(4.5)

\bar{Z}_{n}(\beta):=\frac{1}{M_{x}}\sum_{m=1}^{M_{x}}\frac{\mathrm{det}\big{(}\overline{\mathcal{W}}(\mathbf{x}(\cdot\,;m))\big{)}}{n!\,(2\pi\beta)^{3n/2}\,p(\mathbf{x}(0\,;m))},

which is the Monte Carlo approximation of the multi-dimensional integral

\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\frac{{\rm det}\big{(}\mathcal{\overline{W}}(\mathbf{x})\big{)}}{n!(2\pi\beta)^{3n/2}}\Big{]}{\mathrm{d}}\mathbf{x}_{0}

based on $M_{x}$ samples of independent paths

(4.6)

\mathbf{x}(t\,;i)=\mathbf{B}(t\,;i)+\mathbf{x}(0\,;i),\quad\textup{for }i=1,\dots,M_{x},

where $\mathbf{x}(0\,;i)\in\mathbb{R}^{3n}$ is the $i$ -th sample of the initial position, and $\mathbf{B}(t\,;i)\in\mathbb{R}^{3n}$ is an independent Brownian bridge process sample defined by (2.6), which satisfies $\mathbf{B}(0)=\mathbf{B}(\beta)=0$ . The samples of the initial positions $\{\mathbf{x}(0\,;i)\}_{i=1}^{M_{x}}$ are independently drawn from a mixed normal distribution, with probability density function

(4.7)

p(\mathbf{x})=\frac{1}{2}\Big{(}\frac{1}{\sigma_{1}\sqrt{2\pi}}\,e^{-\frac{|\mathbf{x}|^{2}}{2\sigma_{1}^{2}}}+\frac{1}{\sigma_{2}\sqrt{2\pi}}\,e^{-\frac{|\mathbf{x}|^{2}}{2\sigma_{2}^{2}}}\Big{)}\,,

where $\sigma_{1}^{2}=\beta$ and $\sigma_{2}^{2}=1/\beta$ are parameters adjusting the variance of the sampling distribution for various inverse temperatures $\beta$ . The matrix element in the $k$ -th row and ${\ell}$ -th column of $\mathcal{\overline{W}}(\mathbf{x})$ is given by

	$\displaystyle\mathcal{\overline{W}}_{k\ell}\big{(}\mathbf{x}\big{)}$	$\displaystyle=e^{-\|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)\|^{2}/(2\beta)}e^{-\sum_{m=1}^{M-1}\tilde{V}_{k}(\mathbf{B}_{k}(t_{m})+(1-\frac{t_{m}}{\beta})\mathbf{x}_{k}(0)+\frac{t_{m}}{\beta}\mathbf{x}_{\ell}(0))\Delta t}\times$
		$\displaystyle\quad\times e^{-\big{(}\tilde{V}_{k}(\mathbf{x}_{k}(0))+\tilde{V}_{k}(\mathbf{x}_{\ell}(0))\big{)}\Delta t/2}\,.$

Here we use the trapezoidal method to evaluate the integration in the imaginary time range $[0,\beta]$ numerically, with step size $\Delta t=\beta/M$ , discrete time points $t_{m}=m\,\Delta t$ for $m=0,1,\dots,M$ , and $\mathbf{B}_{k}(t)$ is the $k$ -th component of the Brownian bridge process $\mathbf{B}(t)$ . Due to the absence of pairwise interactions between the fermions in this model, the potential function $V(\mathbf{x},\mathbf{X})$ depends only on the external harmonic oscillator term and is particle-wisely separable

V(\mathbf{x},\mathbf{X})=\sum_{k=1}^{n}\tilde{V}_{k}(\mathbf{x}_{k}),\ \textup{ with }\ \tilde{V}_{k}(\mathbf{x}_{k})=V_{\mathrm{HO}}(\mathbf{x}_{k})=\frac{1}{2}|\mathbf{x}_{k}|^{2}\,.

Specifically, we use different time-step sizes $\Delta t$ to obtain the approximations $\bar{Z}_{n}(\beta)$ for $n=6$ non-interacting fermions at varying $\beta$ values. The results of our numerical implementations are summarized in Table 4.1.

$\beta$	$M_{x}$	$Z_{n}(\beta)$	$\Delta t$	$\bar{Z}_{n}(\beta)$	Relative diff.	Relative CI
$1.0$	$2^{28}$	$1.6978\times 10^{-4}$	$0.025$	$1.6983(8)\times 10^{-4}$	$2.96\,\mathrm{e}(-4)$	$4.54\,\mathrm{e}(-4)$
$1.0$	$2^{28}$	$1.6978\times 10^{-4}$	$0.0125$	$1.6980(8)\times 10^{-4}$	$1.08\,\mathrm{e}(-4)$	$4.56\,\mathrm{e}(-4)$
$1.5$	$2^{28}$	$5.9174\times 10^{-9}$	$0.025$	$5.922(13)\times 10^{-9}$	$8.45\,\mathrm{e}(-4)$	$2.14\,\mathrm{e}(-3)$
$1.5$	$2^{28}$	$5.9174\times 10^{-9}$	$0.0125$	$5.915(13)\times 10^{-9}$	$3.43\,\mathrm{e}(-4)$	$2.14\,\mathrm{e}(-3)$
$2.0$	$2^{28}$	$6.8663\times 10^{-13}$	$0.025$	$7.1(5)\times 10^{-13}$	$3.59\,\mathrm{e}(-2)$	$7.81\,\mathrm{e}(-2)$
$2.0$	$2^{28}$	$6.8663\times 10^{-13}$	$0.0125$	$6.9(5)\times 10^{-13}$	$7.25\,\mathrm{e}(-3)$	$7.65\,\mathrm{e}(-2)$

Table 4.1. Case V1: Approximation of partition function

\bar{Z}_{n}(\beta)

obtained by (4.5) for

n=6

non-interacting fermions in dimension 3, compared with the exact value from (4.2). The digit in the parenthesis of the

\bar{Z}_{n}(\beta)

values denotes the statistical error in our numerical result by the Monte Carlo integration method. The relative difference between

\bar{Z}_{n}(\beta)

and

Z_{n}(\beta)

is computed by

|\bar{Z}_{n}(\beta)-Z_{n}(\beta)|/Z_{n}(\beta)

. The 95% relative confidence interval of

\bar{Z}_{n}(\beta)

is recorded in the last column. The total sample size of the independent Monte Carlo estimators for evaluating the multi-dimensional integral in

\mathbb{R}^{3n}

is denoted by

M_{x}

Furthermore, we also obtain estimations for the system mean-field energy $h(\beta)$ based on the derivatives of the approximated partition function $\bar{Z}_{n}(\beta)$ with respect to the parameter $\beta$

h(\beta)=-\partial_{\beta}\log\big{(}Z_{n}(\beta)\big{)}\simeq-\partial_{\beta}\log{\big{(}\bar{Z}_{n}(\beta)\big{)}}=:\bar{h}(\beta)\,.

We apply a rescaling technique on the sampled Brownian bridge paths to obtain the derivatives of the matrix element $\partial_{\beta}\mathcal{W}_{k\ell}$ . This technique allows a formulation using integration based on standard Brownian bridge samples over the time range $[0,1]$ , avoiding the need for evaluating the derivatives of the Brownian bridge processes with respect to $\beta$ over the time range $[0,\beta]$ . Further details regarding this rescaling technique are provided in the beginning of Section 4.2 and in Appendix A.1.

To establish a reference value for the mean-field energy $h(\beta)$ , we utilize a simple central difference formula to approximate the derivative:

\partial_{\beta}\log\big{(}Z_{n}(\beta)\big{)}\simeq\frac{\log{\big{(}Z_{n}(\beta+\Delta\beta)\big{)}}-\log{\big{(}Z_{n}(\beta-\Delta\beta)\big{)}}}{2\Delta\beta}+\mathcal{O}\big{(}\Delta\beta\,^{2}\big{)}\,.

Here, $Z_{n}(\beta+\Delta\beta)$ and $Z_{n}(\beta-\Delta\beta)$ are computed using the exact recursive formula (4.2), and an inner loop gradually reduces the $\Delta\beta$ value to attain the desired accuracy for the reference value of $h(\beta)$ .

The corresponding results for our approximation of the system mean-field energy $\bar{h}(\beta)$ are summarized in Table 4.2. In Figure 1(b), we plot the error of our Monte Carlo approximations to the partition function $\bar{Z}_{n}(\beta)$ and to the mean-field energy $\bar{h}(\beta)$ with increasing sample size $M_{x}$ at inverse temperature $\beta=1$ .

$\beta$	$M_{x}$	$h(\beta)$	$\Delta t$	$\bar{h}(\beta)$	Relative diff.	Relative CI
$1.0$	$2^{28}$	$22.7799$	$0.025$	$22.778(10)$	$7.81\,\mathrm{e}(-5)$	$4.44\,\mathrm{e}(-4)$
$1.0$	$2^{28}$	$22.7799$	$0.0125$	$22.780(10)$	$2.09\,\mathrm{e}(-5)$	$4.45\,\mathrm{e}(-4)$
$1.5$	$2^{28}$	$18.9572$	$0.025$	$18.95(3)$	$2.14\,\mathrm{e}(-4)$	$1.66\,\mathrm{e}(-3)$
$1.5$	$2^{28}$	$18.9572$	$0.0125$	$18.96(3)$	$4.86\,\mathrm{e}(-5)$	$1.66\,\mathrm{e}(-3)$
$2.0$	$2^{28}$	$17.4894$	$0.025$	$17.2(8)$	$1.53\,\mathrm{e}(-2)$	$4.44\,\mathrm{e}(-2)$
$2.0$	$2^{28}$	$17.4894$	$0.0125$	$17.4(8)$	$5.88\,\mathrm{e}(-3)$	$4.54\,\mathrm{e}(-2)$

Table 4.2. Case V1: Approximation of the mean-field energy

\bar{h}(\beta)

for

n=6

non-interacting fermions in dimension 3, compared with the corresponding exact value

h(\beta)

. The digit in the parenthesis of the

\bar{h}(\beta)

values denotes the statistical error by the Monte Carlo integration method. The relative difference between

\bar{h}(\beta)

and

h(\beta)

is computed by

|\bar{h}(\beta)-h(\beta)|/h(\beta)

, and the 95% relative confidence interval of

\bar{h}(\beta)

is recorded in the last column. The total sample size of the independent Monte Carlo estimators for evaluating the multi-dimensional integral in

\mathbb{R}^{3n}

is denoted by

M_{x}

Refer to caption — (a) Error in $\bar{Z}_{n}(\beta)$ .

From Tables 4.1 and 4.2, we observe that the statistical error of our Monte Carlo estimation of the high-dimensional integration in the Feynman-Kac formulation is dominating the time discretization error since the difference between using $\Delta t=0.025$ and $\Delta t=0.0125$ is smaller than the corresponding confidence interval based on $M_{x}=2^{28}$ samples. The statistical error of the Monte Carlo estimation grows also quickly as the inverse temperature $\beta$ increases, due to the more severe cancellation of permuted particle configurations at higher $\beta$ value.

Moreover, since no additional approximations are introduced for this example with the absence of pairwise interactions between particles, we have that $\bar{Z}_{n}(\beta)$ and $\bar{h}(\beta)$ will approach the exact partition function $Z_{n}(\beta)$ and mean-field energy $h(\beta)$ , in the limit of $\Delta t\to 0$ and $M_{x}\to\infty$ . This is also manifested by the error curves in Figure 1(b), where no biases are observed as the sample size $M_{x}$ increases. Also, the errors decrease with a rate approximately proportional to $M_{x}^{-\frac{1}{2}}$ , confirming the well-known convergence rate of Monte Carlo approximation of the high-dimensional integral. As a comparison, we also calculate the partition function for a system consisting of distinguishable particles. The computationally faster approximation of the fermionic mean-field energy by distinguishable particle partition function gives an error of $15\%$ for $\beta=1$ and $25\%$ for $\beta=1.5$ .

4.2. Electrons with Coulomb interaction under a harmonic trap

For a model including nuclei and involving pairwise Coulomb interactions, our approximation to the system partition function reads

\mathcal{Z}_{\nu}=e^{-\frac{\beta}{2}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}Z_{i}Z_{j}/|\mathbf{X}_{i}-\mathbf{X}_{j}|}\int_{\mathbb{R}^{dn}}\mathbb{E}\big{[}\frac{\det{\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}}}{n!(2\pi\beta)^{dn/2}}\big{]}\,\mathrm{d}\mathbf{x}(0)\,,

as given in (3.9). Recalling the corresponding definition of the matrix element

\mathcal{W}_{k\ell}(\beta,\mathbf{x}):=e^{-|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}/(2\beta)}e^{\int_{0}^{\beta}\sum_{j=1}^{N}Z_{j}/|\mathbf{x}_{k}^{\ell}(t)-\mathbf{X}_{j}|{\mathrm{d}}t}e^{-\frac{1}{2}\int_{0}^{\beta}\sum_{j=1,j\neq k}^{n}1/|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|{\mathrm{d}}t}\,,

where

\nu_{j}=\left\{\begin{array}[]{ll}j,&\mbox{ if }j\neq\ell\,,\\ k,&\mbox{ if }j=\ell\,,\end{array}\right.\,

we can take the derivative of $\mathcal{Z}_{\nu}$ with respect to the inverse temperature $\beta$ , and obtain the corresponding approximation of the system mean-field energy $h$ by

h:=\frac{\mathrm{Tr}\,(H_{e}e^{-\beta H_{e}})}{\mathrm{Tr}\,(e^{-\beta H_{e}})}\approx-\frac{\partial_{\beta}\,\mathcal{Z}_{\nu}(\beta)}{\mathcal{Z}_{\nu}(\beta)}\,.

To evaluate the derivative $\partial_{\beta}\,\mathcal{Z}_{\nu}(\beta)$ , we need to have the derivative of the determinant value $\partial_{\beta}\det{\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}}$ . Jacobi’s formula yields

(4.8)

\frac{\partial}{\partial\beta}\mathrm{det}\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}=\mathrm{Tr}\Big{(}\mathrm{adj}\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}\,\frac{\partial\mathcal{W}(\beta,\mathbf{x})}{\partial\beta}\Big{)}\,,

where $\frac{\partial\mathcal{W}(\beta,\mathbf{x})}{\partial\beta}$ is the element-wise derivative of matrix $\mathcal{W}(\beta,\mathbf{x})$ with respect to $\beta$ , and $\mathrm{adj}\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}$ denotes the adjugate of the matrix $\mathcal{W}(\beta,\mathbf{x})$ , see [29]. The $(i,j)$ -element of this adjugate matrix is defined by the transpose of the co-factor matrix of $\mathcal{W}(\beta,\mathbf{x})$ , i.e.,

\mathrm{adj}\big{(}\mathcal{W}(\beta,\mathbf{x})\big{)}_{ij}:=(-1)^{i+j}\det{\big{(}M_{ji}\big{)}},

with $M_{pq}$ denoting the submatrix of $\mathcal{W}(\beta,\mathbf{x})$ obtained by deleting it $p$ -th row and $q$ -th column.

We consider a system consisting of $n$ identical electrons, with pairwise Coulomb repulsive interactions under a confining harmonic oscillator potential. The matrix element $\mathcal{W}_{k\ell}(\beta,\mathbf{x})$ takes the form

(4.9)

\mathcal{W}_{k\ell}(\beta,\mathbf{x})=e^{-\frac{|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}}{2\beta}}e^{-\int_{0}^{\beta}\tilde{V}_{k\ell}(\mathbf{x}(t)){\mathrm{d}}t}\,,

where

(4.10)

\tilde{V}_{k\ell}\big{(}\mathbf{x}(t)\big{)}=\frac{1}{2}|\mathbf{x}_{k}^{\ell}(t)|^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|}\,,

and the parameter $\lambda$ measures the relative strength of the Coulomb interactions compared to the harmonic oscillator potential. Throughout the numerical experiments of this subsection, we choose the relative force parameter $\lambda=0.5$ . Here $\mathbf{x}_{k}^{\ell}(t)$ denotes a Brownian bridge path with initial position $\mathbf{x}_{k}(0)$ at $t=0$ and ending position $\mathbf{x}_{\ell}(0)$ at $t=\beta$ , defined by

(4.11)

\mathbf{x}_{k}^{\ell}(t):=\mathbf{B}_{k}(t)+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\ell}(0)\,,\ t\in[0,\beta]\,,

where $\mathbf{B}_{k}(t)$ denotes a Brownian bridge with $\mathbf{B}_{k}(0)=\mathbf{B}_{k}(\beta)=0$ , as defined in (2.6).

In order to determine $\frac{\partial\mathcal{W}(\beta,\mathbf{x})}{\partial\beta}$ , it is useful to define the stochastic process $\overline{\mathbf{B}}:[0,1]\times\Omega\to\mathbb{R}^{dn}$ by

(4.12)

\overline{\mathbf{B}}(s):=\frac{1}{\sqrt{\beta}}\mathbf{B}(\beta s)=\frac{1}{\sqrt{\beta}}\mathbf{W}(\beta s)-\frac{s}{\sqrt{\beta}}\mathbf{W}(\beta)\,,

where $\mathbf{W}:[0,\beta]\times\Omega\to\mathbb{R}^{dn}$ is the standard Wiener process in $\mathbb{R}^{dn}$ . The process $\overline{\mathbf{B}}$ defined in (4.12) is a standard Brownian bridge process, as verified in Appendix A.2. Using this re-scaling in the time variable, we have the alternative expression for the Brownian bridge path $\mathbf{x}_{k}^{\ell}(t)$ ,

(4.13)

\mathbf{x}_{k}^{\ell}(t)=\sqrt{\beta}\,\overline{\mathbf{B}}_{k}(\frac{t}{\beta})+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\ell}(0)\,,\ t\in[0,\beta]\,.

To simplify the notation, we define for all indices $k$ and ${\ell}$

	$\displaystyle a_{k\ell}(\beta,\mathbf{x})$	$\displaystyle:=e^{-\frac{\|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)\|^{2}}{2\beta}}\,,$
	$\displaystyle\gamma_{k\ell}(\beta,\mathbf{x})$	$\displaystyle:=\int_{0}^{\beta}\tilde{V}_{k\ell}(\mathbf{x}(t))\,{\mathrm{d}}t\,,$
	$\displaystyle\textup{and\ }\ b_{k\ell}(\beta,\mathbf{x})$	$\displaystyle:=e^{-\gamma_{k\ell}(\beta,\mathbf{x})}\,,$

then the matrix element $\mathcal{W}_{k\ell}(\beta,\mathbf{x})$ in (4.9) and its derivative with respect to $\beta$ can be written as

\mathcal{W}_{k\ell}(\beta,\mathbf{x})=a_{k\ell}(\beta,\mathbf{x})\,b_{k\ell}(\beta,\mathbf{x})\,,

and

(4.14)

\partial_{\beta}\,\mathcal{W}_{k\ell}(\beta,\mathbf{x})=\frac{\partial a_{k\ell}(\beta,\mathbf{x})}{\partial\beta}\,b_{k\ell}(\beta,\mathbf{x})+a_{k\ell}(\beta,\mathbf{x})\,\frac{\partial b_{k\ell}(\beta,\mathbf{x})}{\partial\beta}\,.

Specifically, we have

(4.15)		$\displaystyle\frac{\partial a_{k\ell}(\beta,\mathbf{x})}{\partial\beta}$	$\displaystyle=\exp{\big{(}-\frac{\|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)\|^{2}}{2\beta}\big{)}}\,\frac{\|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)\|^{2}}{2\beta^{2}}\,,$
(4.15)		$\displaystyle\frac{\partial b_{k\ell}(\beta,\mathbf{x})}{\partial\beta}$	$\displaystyle=e^{-\gamma_{k\ell}(\beta,\mathbf{x})}\,\Big{(}-\frac{\partial}{\partial\beta}\,\gamma_{k\ell}(\beta,\mathbf{x})\Big{)}\,.$

Using the standard Brownian bridge process $\overline{\mathbf{B}}(s)$ with time variable $s\in[0,1]$ and introducing the shorthand notation

\bar{\mathbf{x}}_{k}^{\ell}(s):=\sqrt{\beta}\,\overline{\mathbf{B}}_{k}(s)+(1-s)\,\mathbf{x}_{k}(0)+s\,\mathbf{x}_{\ell}(0)\,,\textup{ for }s\in[0,1]\,,

we have the expression

\gamma_{k\ell}(\beta,\mathbf{x})=\int_{0}^{1}\Big{(}\frac{1}{2}\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{|}}\Big{)}\beta\,{\mathrm{d}}s=:\beta\int_{0}^{1}\tilde{V}_{k\ell}\big{(}\bar{\mathbf{x}}(s)\big{)}\mathrm{d}s.

Thus by employing the rescaling technique (4.13), we avoid the derivative $\partial_{\beta}\mathbf{B}(t)$ and obtain

\partial_{\beta}\gamma_{k\ell}(\beta,\mathbf{x})=\int_{0}^{1}\Big{(}\tilde{V}_{k\ell}\big{(}{\bar{\mathbf{x}}(s)}\big{)}+\frac{\sqrt{\beta}}{2}\,\overline{\mathbf{B}}(s)\cdot\nabla\tilde{V}_{k\ell}\big{(}{\bar{\mathbf{x}}(s)}\big{)}\Big{)}\,\mathrm{d}s.

A more detailed expression for the implementation with pairwise Coulomb repulsion and a confining harmonic potential is provided in Appendix A.1.

With the absence of nuclei, the partition function $\mathcal{Z}_{\nu}(\beta)$ becomes

\mathcal{Z}_{\nu}(\beta)=\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\frac{\mathrm{det}(\mathcal{W}(\mathbf{x})}{n!(2\pi\beta)^{3n/2}}\Big{]}\,\mathrm{d}\mathbf{x}\,,

and the derivative of $\mathcal{Z}_{\nu}(\beta)$ with respect to $\beta$ is given by

	$\displaystyle\partial_{\beta}\,\mathcal{Z}_{\nu}(\beta)$	$\displaystyle=\int_{\mathbb{R}^{3n}}\frac{1}{n!}\mathbb{E}\Big{[}\frac{1}{(2\pi\beta)^{3n/2}}\frac{\partial}{\partial\beta}\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}-\frac{3n}{2\beta(2\pi)^{3n/2}}\beta^{-\frac{3n+2}{2}}\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x}$
		$\displaystyle=\int_{\mathbb{R}^{3n}}\frac{1}{n!}\frac{1}{(2\pi\beta)^{3n/2}}\mathbb{E}\Big{[}\mathrm{Tr}\Big{(}\mathrm{adj}\big{(}\mathcal{W}(\mathbf{x})\big{)}\,\frac{\partial\mathcal{W}(\mathbf{x})}{\partial\beta}\Big{)}-\frac{3n}{2\beta}\,\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x}\,.$

Therefore, based on the approximation of the partition function $\mathcal{Z}_{\nu}$ , the system mean-field energy $h$ can be approximated by

(4.16)		$\displaystyle h_{\nu}$	$\displaystyle:=-\frac{\partial_{\beta}\,\mathcal{Z}_{\nu}(\beta)}{\mathcal{Z}_{\nu}(\beta)}$
(4.16)			$\displaystyle=-\frac{\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\mathrm{Tr}\Big{(}\mathrm{adj}\big{(}\mathcal{W}(\mathbf{x})\big{)}\,\frac{\partial\mathcal{W}(\mathbf{x})}{\partial\beta}\Big{)}-\frac{3n}{2\beta}\,\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x_{0}}}{\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x_{0}}}\,,$

where $\frac{\partial\mathcal{W}(\mathbf{x})}{\partial\beta}$ is computed by the element-wise formula (4.14).

In the numerical experiments, the integrals in the $\mathbb{R}^{3n}$ space appearing in (4.16) are estimated with the Monte Carlo integration method by independently sampling initial states $\mathbf{x}(0\,;i)\in\mathbb{R}^{3n}$ for $i=1,2,\cdots,M_{x}$ , and sampling, for each initial state $\mathbf{x}(0\,;i)$ , an associated independent Brownian bridge path $\mathbf{B}(t\,;i)$ in $\mathbb{R}^{3n}$ , which gives a sample for the whole path $\mathbf{x}(t\,;i)$ for $t\in[0,\beta]$ following (4.6). For simplicity of notation, we omit the variable $t$ for the stochastic Brownian paths $\mathbf{x}$ in (4.16). More details on statistical analysis of simulations and estimators of the mean-field approximation $\bar{h}_{v}$ are provided in Appendix B.

Similarly to Section 4.1, our sampling of Brownian bridge paths is based on a discretization scheme in time with step size $\Delta t$ , and for evaluating the element-wise derivative $\frac{\partial\mathcal{W}(\mathbf{x})}{\partial\beta}$ , we use the same numerical integration scheme with trapezoidal rule as shown in (3.4), which gives a discretization error proportional to $\Delta t^{2}$ in our result. Our numerical approximation of mean-field energy involving also the time discretization error is denoted by $\bar{h}_{\nu}$ .

By varying the time-step size $\Delta t$ and taking a large sample size $M_{x}=2^{26}$ for Monte Carlo integral evaluation, we obtain the approximate mean-field energy $\bar{h}_{\nu}$ for a system containing $n=6$ electrons with Coulomb interaction under a confining harmonic trap $V_{\textup{HO}}(\mathbf{x})=\frac{1}{2}|\mathbf{x}|^{2}$ at varying inverse temperature $\beta$ , and make a comparison with the reference result $h_{\textup{ref}}$ obtained by the established algorithm in [11]. The corresponding results are summarized in Table 4.3. As the system inverse temperature increases from $\beta=0.5$ to $\beta=2$ , the relative difference level in our approximation $\bar{h}_{\nu}$ rises from $10^{-4}$ to $10^{-2}$ .

$\beta$	$h_{\textup{ref}}$	$\Delta t$	$\bar{h}_{\nu}$	Relative diff.	Relative CI
$0.5$	$41.66(1)$	$0.025$	$41.655(3)$	$1.14\,\mathrm{e}(-4)$	$6.48\,\mathrm{e}(-5)$
$0.5$	$41.66(1)$	$0.0125$	$41.653(3)$	$1.64\,\mathrm{e}(-4)$	$6.47\,\mathrm{e}(-5)$
$1.0$	$26.692(9)$	$0.025$	$26.711(4)$	$7.12\,\mathrm{e}(-4)$	$1.48\,\mathrm{e}(-4)$
$1.0$	$26.692(9)$	$0.0125$	$26.716(4)$	$8.87\,\mathrm{e}(-4)$	$1.38\,\mathrm{e}(-4)$
$1.5$	$22.63(7)$	$0.025$	$22.774(8)$	$6.37\,\mathrm{e}(-3)$	$3.34\,\mathrm{e}(-4)$
$1.5$	$22.63(7)$	$0.0125$	$22.774(7)$	$6.37\,\mathrm{e}(-3)$	$3.23\,\mathrm{e}(-4)$
$2.0$	$22.1(5)$	$0.025$	$20.66(5)$	$6.52\,\mathrm{e}(-2)$	$2.47\,\mathrm{e}(-3)$
$2.0$	$22.1(5)$	$0.0125$	$20.69(6)$	$6.39\,\mathrm{e}(-2)$	$2.69\,\mathrm{e}(-3)$

Table 4.3. Case V2: The approximate mean-field energy

\bar{h}_{\nu}

obtained by Monte Carlo approximation of (4.16) with sample size

M_{x}=2^{26}

for

d=3

n=6

, compared to the reference value from [11]. The digit in the parenthesis of

\bar{h}_{\nu}

values denotes the statistical error from the Monte Carlo integration method. The relative difference between

\bar{h}_{\nu}

and

h_{\mathrm{ref}}

is computed by

|\bar{h}_{\nu}-h_{\mathrm{ref}}|/h_{\mathrm{ref}}

, and the last column records the 95% relative confidence interval of

\bar{h}_{\nu}

Moreover, we implement our path integral method with Feymann-Kac formulation on a $2$ -dimensional quantum dot model, considering pair-wise Coulomb repulsion between electrons and still the confining harmonic trap $V_{\mathrm{HO}}$ . The system inverse temperature is fixed to be $\beta=0.3$ and $\beta=1$ , respectively, while the number of electrons increases from $n=6$ to $20$ and from $n=3$ to $10$ . This $2$ -dimensional model is also studied in [11] with the reference values for the mean-field energy $h_{\mathrm{ref}}$ for benchmarking our approximation $\bar{h}_{\nu}$ . We summarize our numerical results obtained with different time step sizes $\Delta t$ in Table 4.4.

$\beta$	$n$	$h_{\textup{ref}}$	$\Delta t$	$\bar{h}_{\nu}$	Relative diff.	Relative CI
$1$	$3$	$8.719(3)$	$0.025$	$8.717(3)$	$2.74\,\mathrm{e}(-4)$	$3.83\,\mathrm{e}(-4)$
	$3$	$8.719(3)$	$0.0125$	$8.717(3)$	$2.75\,\mathrm{e}(-4)$	$3.95\,\mathrm{e}(-4)$
	$6$	$22.82(5)$	$0.025$	$22.82(3)$	$2.13\,\mathrm{e}(-4)$	$1.20\,\mathrm{e}(-3)$
	$6$	$22.82(5)$	$0.0125$	$22.79(2)$	$1.03\,\mathrm{e}(-3)$	$1.04\,\mathrm{e}(-3)$
	$10$	$49(3)$	$0.025$	$48.3(5)$	$1.47\,\mathrm{e}(-2)$	$9.92\,\mathrm{e}(-3)$
	$10$	$49(3)$	$0.0125$	$48.5(4)$	$1.04\,\mathrm{e}(-2)$	$8.89\,\mathrm{e}(-3)$
$0.3$	$6$	$46.45(1)$	$0.025$	$46.44(1)$	$2.17\,\mathrm{e}(-4)$	$2.83\,\mathrm{e}(-4)$
	$6$	$46.45(1)$	$0.0125$	$46.45(1)$	$9.79\,\mathrm{e}(-6)$	$2.85\,\mathrm{e}(-4)$
	$10$	$84.92(4)$	$0.025$	$84.90(2)$	$2.69\,\mathrm{e}(-4)$	$2.85\,\mathrm{e}(-4)$
	$10$	$84.92(4)$	$0.0125$	$84.89(2)$	$3.04\,\mathrm{e}(-4)$	$2.86\,\mathrm{e}(-4)$
	$20$	$203(1)$	$0.025$	$203.5(2)$	$2.63\,\mathrm{e}(-3)$	$1.04\,\mathrm{e}(-3)$
	$20$	$203(1)$	$0.0125$	$203.4(2)$	$2.07\,\mathrm{e}(-3)$	$8.85\,\mathrm{e}(-4)$

Table 4.4. Case V2: The approximate mean-field energy

\bar{h}_{\nu}

obtained with equation (4.16) using sample size

M_{x}=2^{22}

for

d=2

, compared with the reference value from [11]. The digit in the parenthesis of

\bar{h}_{\nu}

data denotes the statistical error from the Monte Carlo integration method. The relative difference between

\bar{h}_{\nu}

and

h_{\mathrm{ref}}

is computed by

|\bar{h}_{\nu}-h_{\mathrm{ref}}|/h_{\mathrm{ref}}

, and the last column records the 95% relative confidence interval of

\bar{h}_{\nu}

For both test cases with $\beta=1$ and $\beta=0.3$ , our approximation $\bar{h}_{\nu}$ employs $M_{x}=2^{22}$ samples for Monte Carlo integral evaluation and agrees with the reference value $h_{\mathrm{ref}}$ on their corresponding level of statistical confidence intervals.

To further survey the bias level of our mean-field energy approximation $\bar{h}_{\nu}$ in the test case V2, we compare our results with both the reference value $h_{\mathrm{ref}}$ in [11], and with a Monte Carlo approximation $h_{\mathrm{tensor}}$ of the exact mean-field energy based on the tensor formulation (3.6), for $n=6$ fermions in dimension $d=3$ . The sample size $M_{x}$ employed for the Monte Carlo integral evaluation of $\bar{h}_{\nu}$ is $2^{26}$ . For the Monte Carlo integral approximation of the more computationally demanding tensor formula, however, the sample size is taken to be $2^{22}$ , which leads to relatively large confidence intervals for the obtained $h_{\mathrm{tensor}}$ data. The numerical results for $h_{\mathrm{ref}}$ , $h_{\mathrm{tensor}}$ , and $\bar{h}_{\nu}$ at different $\beta$ values are summarized in Table 4.5, from which we observe that the relative difference $|\bar{h}_{\nu}-h_{\mathrm{tensor}}|/h_{\mathrm{tensor}}$ increases as $\beta$ increases. Also, the differences $|\bar{h}_{\nu}-h_{\mathrm{tensor}}|$ and $|h_{\mathrm{ref}}-h_{\mathrm{tensor}}|$ are smaller than the statistical uncertainty of $h_{\mathrm{tensor}}$ , suggesting that the statistical error from the high-dimensional Monte Carlo integral evaluation is still dominating.

$\beta$	$h_{\textup{ref}}$	$\Delta t$	$h_{\mathrm{tensor}}$	$\bar{h}_{\nu}$	$\frac{\|h_{\mathrm{tensor}}-\bar{h}_{\nu}\|}{h_{\mathrm{tensor}}}$	Relative CI of $h_{\mathrm{tensor}}$
$0.5$	$41.66(1)$	$0.025$	$41.65(7)$	$41.655(3)$	$7.20\,\mathrm{e}(-5)$	$1.64\,\mathrm{e}(-3)$
$0.5$	$41.66(1)$	$0.0125$	$41.66(7)$	$41.653(3)$	$3.87\,\mathrm{e}(-5)$	$1.64\,\mathrm{e}(-3)$
$1.0$	$26.692(9)$	$0.025$	$26.7(1)$	$26.711(4)$	$4.12\,\mathrm{e}(-4)$	$5.71\,\mathrm{e}(-3)$
$1.0$	$26.692(9)$	$0.0125$	$26.7(1)$	$26.716(4)$	$5.99\,\mathrm{e}(-4)$	$5.24\,\mathrm{e}(-3)$
$1.5$	$22.63(7)$	$0.025$	$22.8(2)$	$22.774(8)$	$1.13\,\mathrm{e}(-3)$	$1.05\,\mathrm{e}(-2)$
$1.5$	$22.63(7)$	$0.0125$	$22.8(2)$	$22.774(7)$	$1.09\,\mathrm{e}(-3)$	$8.59\,\mathrm{e}(-3)$

Table 4.5. Case V2: The approximate mean-field energy

\bar{h}_{\nu}

obtained by Monte Carlo approximation of (4.16) with sample size

M_{x}=2^{26}

for

d=3

n=6

, compared to the reference value from [11] and from the tensor formula (3.6) with sample size

2^{22}

. The digit in the parenthesis of

h_{\mathrm{tensor}}

and

\bar{h}_{\nu}

values denotes the statistical error from the Monte Carlo integration method. The relative difference between

\bar{h}_{\nu}

and

h_{\mathrm{tensor}}

is computed by

|\bar{h}_{\nu}-h_{\mathrm{tensor}}|/h_{\mathrm{tensor}}

, and the last column records the 95% relative confidence interval of

h_{\mathrm{tensor}}

In addition, we implement the perturbation of the determinant formulation following (3.14), (3.15) and (3.16), with the aim of obtaining a sensitivity study for the $\bar{h}_{\nu}$ approximation. We compute first the approximated partition function $\bar{Z}_{\mathrm{perturb}}(\beta)$ , using the expected values of $\mathrm{det}\big{(}\widetilde{\mathcal{W}}\big{)}$ based on 100 independent samples of perturbed Brownian paths, and then apply a central difference quotient formula to approximate the mean-field value $\bar{h}_{\mathrm{perturb}}(\beta)$ as

(4.17)

\bar{h}_{\mathrm{perturb}}(\beta):=\frac{\bar{Z}_{\mathrm{perturb}}(\beta+\Delta\beta)-\bar{Z}_{\mathrm{perturb}}(\beta-\Delta\beta)}{2\Delta\beta}.

The results of $\bar{h}_{\mathrm{perturb}}$ and the non-perturbed approximation $\bar{h}_{\nu}$ using the determinant formulation are summarized in Table 4.6. These results are then compared with the reference mean-field value $h_{\mathrm{ref}}$ . Specifically, for $n=3$ , we obtain $h_{\mathrm{ref}}$ using the tensor formula (3.6), while for $n=6$ , we rely on the result from [11] to obtain a reliable reference value $h_{\mathrm{ref}}$ . The relative error of the perturbed mean-field value $\bar{h}_{\mathrm{perturb}}$ is at most one order of magnitude higher than the relative error of $\bar{h}_{\nu}$ derived from the determinant formulation (4.16), which indicates that $\bar{h}_{\mathrm{perturb}}$ can be used to numerically roughly estimate the accuracy of $\bar{h}_{\nu}$ without using a reference value. In both test cases with $\beta=1.0$ and $\beta=1.5$ , we employ $M_{x}=2^{22}$ samples for the Monte Carlo evaluation of the integral in the high-dimensional $\mathbb{R}^{3n}$ space, and for the central difference formula (4.17) we use $\Delta\beta=0.01$ .

$\beta$	$n$	$h_{\textup{ref}}$	$\Delta t$	$\bar{h}_{\nu}$	$\frac{\|h_{\mathrm{ref}}-\bar{h}_{\nu}\|}{h_{\mathrm{ref}}}$	$\bar{h}_{\mathrm{perturb}}$	$\frac{\|h_{\mathrm{ref}}-\bar{h}_{\mathrm{perturb}}\|}{h_{\mathrm{ref}}}$
$1.0$	$3$	$11.355(3)$	$0.025$	$11.356(3)$	$1.25\,\mathrm{e}(-4)$	$11.337(4)$	$1.55\,\mathrm{e}(-3)$
	$3$	$11.355(3)$	$0.0125$	$11.356(3)$	$7.65\,\mathrm{e}(-5)$	$11.337(3)$	$1.63\,\mathrm{e}(-3)$
	$6$	$26.692(9)$	$0.025$	$26.72(1)$	$1.52\,\mathrm{e}(-4)$	$26.60(1)$	$3.40\,\mathrm{e}(-3)$
	$6$	$26.692(9)$	$0.0125$	$26.71(1)$	$1.08\,\mathrm{e}(-4)$	$26.60(1)$	$3.56\,\mathrm{e}(-3)$
$1.5$	$3$	$9.157(2)$	$0.025$	$9.163(4)$	$6.23\,\mathrm{e}(-4)$	$9.129(5)$	$3.02\,\mathrm{e}(-3)$
	$3$	$9.157(2)$	$0.0125$	$9.159(4)$	$1.86\,\mathrm{e}(-4)$	$9.132(5)$	$2.71\,\mathrm{e}(-3)$
	$6$	$22.63(7)$	$0.025$	$22.83(3)$	$8.79\,\mathrm{e}(-3)$	$22.60(5)$	$1.20\,\mathrm{e}(-3)$
	$6$	$22.63(7)$	$0.0125$	$22.85(2)$	$9.87\,\mathrm{e}(-3)$	$22.55(5)$	$3.55\,\mathrm{e}(-3)$

Table 4.6. Case V2: The approximate mean-field energy

\bar{h}_{\nu}

for

d=3

, obtained with equation (4.16) and the corresponding result

\bar{h}_{\mathrm{perturb}}

obtained with the perturbation determinant formula (4.17) using sample size

M_{x}=2^{22}

, compared with the reference value

h_{\mathrm{ref}}

obtained by the tensor formula based on (3.6) and from the benchmark result in [11]. The digit in the parenthesis of

\bar{h}_{\nu}

and

\bar{h}_{\mathrm{perturb}}

denotes the statistical error from the Monte Carlo integration method. All the values in this table are computed with parameter

c_{*}=2

, which is determined according to the cycle length for the given number of particles

n

. Numerical tests with

c_{\ast}

varying from 0.6 to 2 indicate an insensitivity with respect to this parameter, since

|\bar{h}_{\mathrm{perturb}}-h_{\mathrm{ref}}|/h_{\mathrm{ref}}

remains smaller than

5\,\mathrm{e}(-3)

, which is within the range of the statistical error.

5. Conclusions

We develop and analyze a new computational method for approximating the mean-field Hamiltonian of molecular dynamics suitable for estimation of canonical quantum observables in nuclei-electrons systems. The method is based on approximation of path integrals for fermions using Monte Carlo estimation of particular determinants on Brownian bridge paths. In the presented numerical benchmarks we observe that the derived determinant formulation is surprisingly accurate. More specifically, the reported test examples show an approximation error at the same level as the statistical error in the reference solutions taken from [11]. The approximation of mean-feld value $\bar{h}_{\nu}$ , based on discretizations of (4.16), is also consistent with the numerical results of $\bar{h}_{\mathrm{perturb}}$ , obtained from the sensitivity analysis (3.14), (3.15), (3.16), (4.17), shown in Table 4.6. The perturbed mean-field value $\bar{h}_{\mathrm{perturb}}$ thereby provides a computational error indicator for $\bar{h}_{\nu}$ when no reference solution is available. It remains to explain theoretically the accuracy of the approximation based on the formulation (3.9) from the first principles. We remark that similar to other path dependent methods based on the Feynman-Kac formulation the presented method exhibits an increasing bias as the temperature decreases.

For the test case V1 with non-interacting Fermi gas under an external harmonic oscillator potential $V_{\mathrm{HO}}(\mathbf{x})=\frac{1}{2}|\mathbf{x}|^{2}$ , we compare our numerical approximation of the partition function $\bar{Z}_{n}(\beta)$ and the mean-field energy $\bar{h}(\beta)$ with their corresponding exact values. Since the potential $V_{\mathrm{HO}}$ is particle-wise separable, our formulation (3.4) is asymptotically exact as the Monte Carlo integration sample size $M_{x}\to\infty$ and the time step size $\Delta t\to 0$ . From Tables 4.1, 4.2, and Figure 1(b), we observe no biases in the relative error of the approximations as $M_{x}$ increases, and the statistical error from the Monte Carlo approximation of the high-dimensional integral in the Feynman-Kac formulation (3.4) is dominating. Meanwhile, the time discretization error related to $\Delta t$ is smaller than the statistical uncertainty.

For the test case V2, which involves the pairwise Coulomb interactions between electrons and a confining harmonic trap, our numerical approximation to the mean-field energy $\bar{h}_{\nu}$ achieves comparable accuracy to the reference result in [11] for both the two models with dimensions $d=3$ and $d=2$ , as shown in the Tables 4.3 and 4.4. Specifically, as the inverse temperature $\beta$ and the number of fermions $n$ increases, the level of the statistical uncertainty for both methods rises. From Table 4.3 it is observed that the relative difference to the reference mean-field $|\bar{h}_{\nu}-h_{\mathrm{ref}}|/h_{\mathrm{ref}}$ increases as $\beta$ increases, indicating an elevated bias level in the formulation (3.9) for a lower system temperature.

In comparison with the reference result $h_{\mathrm{tensor}}$ using the permutation formula (3.6) we still observe from Table 4.5 that the statistical error from the Monte Carlo integration method is dominating, while for a system comprising $n$ electrons, the computational workload of the tensor formula $h_{\mathrm{tensor}}$ is $\mathcal{O}(n!)$ , however, our approximation $\bar{h}_{\nu}$ only requires a computational complexity proportional to $n^{3}$ .

Acknowledgment

This research was supported by Swedish Research Council grant 2019-03725 and KAUST grant OSR-2019-CRG8-4033.3. The work of P.P. was supported in part by the U.S. Army Research Office Grant W911NF-19-1-0243. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at PDC Center for High Performance Computing, KTH Royal Institute of Technology, partially funded by the Swedish Research Council through grant agreement no. 2022-06725.

References

[1] Anderson J. B., A random-walk simulation of the Schrödinger equation: $H_{3}^{+}$ , J. Chem. Phys. 63, 1499-1503 (1975).
[2] Blankenbecler R., Scalapino D. J. and Sugar R. L., Monte Carlo calculations of coupled boson-fermion systems. I Phys. Rev. D 24, 2278 (1981).
[3] Boninsegni M., Prokof’ev N.V., Svistunov B.V., Worm algorithm and diagrammatic Monte Carlo: a new approach to continuous-space path integral Monte Carlo simulations, Phys. Rev. E Stat. Nonlin. Soft. Matter Phys. 74 036701 (2006).
[4] Borrmann P., and Franke G., Recursion formulas for quantum statistical partition functions, J. Chem. Phys. 98, 2484 (1993).
[5] Calcavecchia F., and Holzmann M., Fermion sign problem in imaginary-time projection continuum quantum Monte Carlo with local interaction, Physical review. E, 93, 043321 (2016).
[6] Cancés E., Defranceschi M., Kutzelnigg W., Le Bris C. and Maday Y., Computational quantum chemistry: a primer, in: Handbook of numerical analysis. Volume X: computational chemistry, Ph. Ciarlet and C. Le Bris eds. (North-Holland, 2003) 3-270.
[7] Cancés E., Jourdain B., and LeLiévre T., Quantum Monte Carlo simulations of fermions: a mathematical analysis of the fixed-node approximation, Mathematical Models and Methods in Applied Sciences, 16, 09, 1403-1440 (2006).
[8] Ceperley D.M., An Overview of Quantum Monte Carlo Methods, Reviews in Mineralogy and Geochemistry, 71 (1): 129-135 (2010).
[9] Ceperley D. M., Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279–355 (1995).
[10] Ceperley, D.M. Fermion nodes. J Stat Phys 63, 1237-1267 (1991).
[11] Dornheim T.,Fermion sign problem in path integral Monte Carlo simulations: Quantum dots, ultracold atoms, and warm dense matter, Phys. Rev. E, Aug 2:(100):023307 (2019).
[12] Dornheim T., Groth S., Filinov A.V. and Bonitz M., Path integral Monte Carlo simulation of degenerate electrons: Permutation-cycle properties, J. Chem. Phys. 151 014108 (2019).
[13] Dumas W. M., On Conditional Wiener Integrals and a Novel Approach to the Fermion Sign Problem, Ph.D thesis at the University of Leicester (2010).
[14] Dumas W.M. and Tretyakov M.V.,Computing conditional Wiener integrals of functionals of a general form, IMA Journal of Numerical Analysis, 31, 1217-1251 (2011).
[15] Feynman R.F. and Hibbs A.R., Quantum Mechanics and Path Integrals, Dover Publications (2010), emended republication of the original published 1965.
[16] Gut A., Probability: A Graduate Course, second ed., Springer New York, NY (2013).
[17] Habershon S., Manolopoulos D.E., Markland T.E. and Miller T.F. 3rd., Ring-polymer molecular dynamics: quantum effects in chemical dynamics from classical trajectories in an extended phase space. Annu Rev Phys Chem., 64 (2013) 387-413.
[18] Hillar C. J. and Lim L.H., Most tensor problems are NP-hard, Journal of the ACM, 60(6), (2013) 1-39.
[19] Hirshberg B., Invernizzi M. and Parrinello M.,Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality, J Chem Phys. May 7;152(17):171102. (2020)
[20] Huang, X., Plechac P., Sandberg M. and Szepessy A, Canonical mean-field molecular dynamics derived from quantum mechanics, ESAIM Math. Model. Numer. Anal. 56 (2022), no. 6, 2197-2238..
[21] Karatzas I. and Shreve S.E., Brownian Motion and Stochastic Calculus, Springer Verlag (1998).
[22] LeBris C., Computational chemistry from the perspective of numerical analysis, Acta Numerica, vol. 14, pp. 363–444, CUP, 2005.
[23] Li Z.-X. and Yao H., Sign-Problem-Free Fermionic Quantum Monte Carlo:Developments and Applications, Annual Review of Condensed Matter Physics 10, 337 (2019)
[24] Lieb E. and Seiringer R., The Stability of Matter in Quantum Mechanics, Cambridge University Press (2010).
[25] Lyubartsev A.P., Simulation of excited states and the sign problem in the path integral Monte Carlo method, J. Phys. A: Math. Gen. 38 (2005) 6659-6674.
[26] Marx D. and Hutter J., Ab Initio Molecular Dynamics: Basic theory and advanced methods, Cambridge University Press (2009).
[27] Perez A., Tuckerman M.E. and Müser M.H., A comparative study of the centroid and ring-polymer molecular dynamics methods for approximating quantum time correlation functions from path integrals, J. Chem. Phys. 2009 May 14;130(18):184105.
[28] Schmidt H.-J., and Schnack J., Thermodynamic fermion-boson symmetry in harmonic oscillator potentials, Physica A-statistical Mechanics and Its Applications 265 (1998): 584-589.
[29] Stewart G.W., On the adjugate matrix, Linear Algebra and its Applications 283 (1998): 151-164.
[30] Shumway J. and Ceperley D. M., Path integral Monte Carlo simulations for fermion systems: Pairing in the electron-hole plasma, J. Phys. IV France 10 (2000) Pr5-3-Pr5-16.
[31] Takahashi M. and Imada M.,Monte Carlo calculation of quantum systems, J. Phys. Soc. Jpn. 53, pp. 963-974 (1984).
[32] Zwanzig R., Nonequilibrium Statistical Mechanics, New York Oxford Univ Press (2001).
[33] Zworski M., Semiclassical Analysis, Providence, RI, American Mathematical Society (2012).
[34] https://github.com/XinHuang2022/Path_Integral_MD_mean_field.git

Appendix A Supplementary for Numerical Implementation

A.1. The derivative of the matrix element $\mathcal{W}_{k\ell}(\beta,\mathbf{x})$ with respect to $\beta$

For a system comprising indistinguishable fermions with Coulomb repulsive interactions and under a confining harmonic oscillator potential, we have the formula for the matrix element $\mathcal{W}_{k\ell}(\beta,\mathbf{x})$ , based on (4.9) and (4.10):

\mathcal{W}_{k\ell}(\beta,\mathbf{x})=e^{-\frac{|\mathbf{x}_{k}(0)-\mathbf{x}_{\ell}(0)|^{2}}{2\beta}}e^{-\int_{0}^{\beta}\tilde{V}_{k,\ell}(\mathbf{x}(t)){\mathrm{d}}t}\,,

where

\tilde{V}_{k,\ell}\big{(}\mathbf{x}(t)\big{)}=\frac{1}{2}|\mathbf{x}_{k}^{\ell}(t)|^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|},

Using the scaling formula (4.13) we can first rewrite the expression for $\gamma_{k,\ell}(\beta,\mathbf{x})$ with a change of the integration variable as

(A.1)

\begin{split}\gamma_{k,\ell}(\beta,\mathbf{x})&=\int_{0}^{\beta}\tilde{V}_{k,\ell}(\mathbf{x}(t))\,{\mathrm{d}}t\\ &=\int_{0}^{\beta}\frac{1}{2}|\mathbf{x}_{k}^{\ell}(t)|^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2|\mathbf{x}_{k}^{\ell}(t)-\mathbf{x}_{j}^{\nu_{j}}(t)|}\,\mathrm{d}t\\ &\dots\big{\{}\ \textup{Plugging in the scaling formulae \eqref{scaling_BB_kell}\,}\ \big{\}}\\ &=\int_{0}^{\beta}\frac{1}{2}\,\big{|}\underbrace{\sqrt{\beta}\,\bar{\mathbf{B}}_{k}(\frac{t}{\beta})+(1-\frac{t}{\beta})\mathbf{x}_{k}(0)+\frac{t}{\beta}\mathbf{x}_{\ell}(0)}_{=:\bar{\mathbf{x}}_{k}^{\ell}(\frac{t}{\beta})}\big{|}^{2}\,\mathrm{d}t\\ \qquad\;\;\;&+\int_{0}^{\beta}\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\frac{\mathrm{d}t}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(\frac{t}{\beta})-\bar{\mathbf{x}}_{j}^{\nu_{j}}(\frac{t}{\beta})\big{|}}\,\\ &\dots\big{\{}\ \textup{Taking the change of variable }\,s:=\frac{t}{\beta},\ \mathrm{d}t=\beta\,\mathrm{d}s\ \big{\}}\\ &=\quad\int_{0}^{1}\frac{1}{2}\,\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{|}^{2}\beta\,\mathrm{d}s+\int_{0}^{1}\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\frac{1}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{|}}\,\beta\,\mathrm{d}s.\\ \end{split}

We have

(A.2)

\frac{\partial}{\partial\beta}\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{|}^{2}=2\,\bar{\mathbf{x}}_{k}^{\ell}(s)\cdot\frac{\partial}{\partial\beta}\bar{\mathbf{x}}_{k}^{\ell}(s)=\bar{\mathbf{x}}_{k}^{\ell}(s)\cdot\frac{1}{\sqrt{\beta}}\bar{\mathbf{B}}_{k}(s)\,,

and

(A.3)

\begin{split}&\frac{\partial}{\partial\beta}\frac{1}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{|}}\\ &=-\frac{1}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{|}^{3}}\,\big{(}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{)}\cdot\frac{\partial}{\partial\beta}\big{(}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{)}\\ &=-\frac{1}{\big{|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{|}^{3}}\,\big{(}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{)}\cdot\Big{(}\frac{1}{2\sqrt{\beta}}\bar{\mathbf{B}}_{k}(s)-\frac{1}{2\sqrt{\beta}}\bar{\mathbf{B}}_{j}(s)\Big{)}\,.\end{split}

Based on equations (A.1), (A.2), and (A.3), we further obtain

(A.4)			$\displaystyle\frac{\partial}{\partial\beta}\,\gamma_{k,\ell}(\beta,\mathbf{x})=\frac{\partial}{\partial\beta}\Big{(}\int_{0}^{1}\big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\beta\,{\mathrm{d}}s\Big{)}$
			$\displaystyle\begin{aligned} =&\int_{0}^{1}\big{(}\frac{1}{2}\,\frac{\partial}{\partial\beta}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\partial}{\partial\beta}\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\beta\,{\mathrm{d}}s\\ &+\int_{0}^{1}\big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\,{\mathrm{d}}s\\ \end{aligned}$
			$\displaystyle\begin{aligned} =&\int_{0}^{1}\Big{(}\frac{1}{2}\,\bar{\mathbf{x}}_{k}^{\ell}(s)\cdot\frac{1}{\sqrt{\beta}}\bar{\mathbf{B}}_{k}(s)\\ &-\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\big{(}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{)}}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}^{3}}\cdot\frac{1}{2\sqrt{\beta}}\big{(}\bar{\mathbf{B}}_{k}(s)-\bar{\mathbf{B}}_{j}(s)\big{)}\Big{)}\beta\,{\mathrm{d}}s\\ &+\int_{0}^{1}\Big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\Big{)}\,{\mathrm{d}}s\\ \end{aligned}$
			$\displaystyle\begin{aligned} =&\int_{0}^{\beta}\Big{(}\frac{1}{2}\,{\mathbf{x}}_{k}^{\ell}(t)\cdot\frac{1}{\beta}{\mathbf{B}}_{k}(t)\\ &-\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\big{(}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{)}}{\big{\|}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{\|}^{3}}\cdot\frac{1}{2\beta}\big{(}{\mathbf{B}}_{k}(t)-{\mathbf{B}}_{j}(t)\big{)}\Big{)}\,{\mathrm{d}}t\\ &+\int_{0}^{\beta}\Big{(}\frac{1}{2}\big{\|}{\mathbf{x}}_{k}^{\ell}(t)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{\|}}\Big{)}\frac{1}{\beta}\,{\mathrm{d}}t\,.\\ \end{aligned}$

Now combining equations (A.4), (4.14), and (4.15), we obtain at the derivative of the matrix element $\partial_{\beta}\,\mathcal{W}_{k,\ell}(\beta,\mathbf{x})$ .

A.2. The time-rescaling technique on the Brownian bridge process $\mathbf{B}(t)$

Applying the covariance function of the standard Wiener process $\mathbf{W}:\mathbb{R}^{+}\times\Omega\to\mathbb{R}$ ,

\mathrm{Cov}\big{(}\mathbf{W}(\tau_{1}),\mathbf{W}(\tau_{2})\big{)}=\min{(\tau_{1},\tau_{2})}\,,

we obtain the covariance function of $\overline{\mathbf{B}}(s)$ as defined in (4.12) for two time values $\tau_{1},\tau_{2}\in[0,1]$ as

		$\displaystyle\mathrm{Cov}\big{(}\overline{\mathbf{B}}(\tau_{1}),\overline{\mathbf{B}}(\tau_{2})\big{)}$
		$\displaystyle=\mathrm{Cov}\big{(}\frac{1}{\sqrt{\beta}}\mathbf{W}(\tau_{1}\beta)-\frac{\tau_{1}}{\sqrt{\beta}}\mathbf{W}(\beta),\frac{1}{\sqrt{\beta}}\mathbf{W}(\tau_{2}\beta)-\frac{\tau_{2}}{\sqrt{\beta}}\mathbf{W}(\beta)\big{)}$
		$\displaystyle=\frac{1}{\beta}\Big{(}\mathrm{Cov}\big{(}\mathbf{W}(\tau_{1}\beta),\mathbf{W}(\tau_{2}\beta)\big{)}-\tau_{1}\mathrm{Cov}\big{(}\mathbf{W}(\tau_{2}\beta),\mathbf{W}(\beta)\big{)}-\tau_{2}\mathrm{Cov}\big{(}\mathbf{W}(\tau_{1}\beta),\mathbf{W}(\beta)\big{)}\Big{)}$
		$\displaystyle\quad+\frac{1}{\beta}\,\tau_{1}\tau_{2}\mathrm{Cov}\big{(}\mathbf{W}(\beta),\mathbf{W}(\beta)\big{)}$
		$\displaystyle=\frac{1}{\beta}\Big{(}\min{(\tau_{1}\beta,\tau_{2}\beta)}-\tau_{1}\tau_{2}\beta-\tau_{2}\tau_{1}\beta+\tau_{1}\tau_{2}\beta\Big{)}$
		$\displaystyle=\min{(\tau_{1},\tau_{2})}-\tau_{1}\tau_{2}\,,$

which satisfies the definition of a standard Brownian bridge process defined with the time range $s\in[0,1]$ . The Brownian bridge $\mathbf{B}_{k}(t)$ in the equation (4.11) can thus be expressed with a scaling factor $\beta$ as

(A.5)

\mathbf{B}_{k}(t)\overset{\mathrm{d}}{=}\sqrt{\beta}\,\overline{\mathbf{B}}_{k}(\frac{t}{\beta})\,,\ t\in[0,\beta]\,,

where $\overline{\mathbf{B}}_{k}(s)$ denotes a standard Brownian bridge starting and ending at the same point $\mathbf{x}_{k}(0)$ with the time variable $s\in[0,1]$ , and ‘ $\mathrm{d}$ ’ denotes equality of the two random processes in distribution.

Appendix B Statistical confidence interval of the approximation of mean-field $\bar{h}_{\nu}$

Based on the approximation (3.9) of the partition function $\mathcal{Z}_{\nu}$ , we have for the mean-field energy the approximate formula (4.16)

h_{\nu}=-\frac{\partial_{\beta}\,\mathcal{Z}_{\nu}(\beta)}{\mathcal{Z}_{\nu}(\beta)}=-\frac{\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\mathrm{Tr}\Big{(}\mathrm{adj}\big{(}\mathcal{W}(\mathbf{x})\big{)}\,\frac{\partial\mathcal{W}(\mathbf{x})}{\partial\beta}\Big{)}-\frac{3n}{2\beta}\,\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x_{0}}}{\int_{\mathbb{R}^{3n}}\mathbb{E}\Big{[}\mathrm{det}\big{(}\mathcal{W}(\mathbf{x})\big{)}\Big{]}\,\mathrm{d}\mathbf{x_{0}}}=:\frac{\mathbb{E}[A]}{\mathbb{E}[B]}\,,

where $A$ and $B$ are the random variables corresponding to the Monte Carlo approximation of the integration in the numerator and the denominator, respectively. For evaluating the Monte Carlo integrals, a sample set of size $N$ is utilized, which yields the approximation of the mean-filed $\bar{h}_{\nu,N}$

(B.1)

\bar{h}_{\nu,N}=\frac{\frac{1}{N}\sum_{n=1}^{N}A_{n}}{\frac{1}{N}\sum_{n=1}^{N}B_{n}}\,,

where $\{A_{n}\}_{n=1}^{N}$ and $\{B_{n}\}_{n=1}^{N}$ are sample sets based on the identically and independently drawn samples of paths $\{\mathbf{x}^{(n)}\}_{n=1}^{N}$ with initial position $\mathbf{x}^{(n)}_{0}$ following the mixed normal distribution with probability density $p(x)$ , which is described by (4.7). Specifically we have

	$\displaystyle A_{n}$	$\displaystyle=\Big{(}\mathrm{Tr}\big{(}\mathrm{adj}\big{(}\mathcal{W}(\mathbf{x}^{(n)})\big{)}\,\frac{\partial\mathcal{W}(\mathbf{x}^{(n)})}{\partial\beta}\big{)}-\frac{3n}{2\beta}\,\mathrm{det}\big{(}\mathcal{W}(\mathbf{x}^{(n)})\big{)}\Big{)}/p(\mathbf{x}^{(n)}_{0})\,,$
	$\displaystyle B_{n}$	$\displaystyle=\mathrm{det}\big{(}\mathcal{W}(\mathbf{x}^{(n)})\big{)}/p(\mathbf{x}^{(n)}_{0})\,.$

In this section, we derive estimates of the above statistical approximation of $h_{\nu}$ based on the Berry-Esseen theorem.

Applying the central limit theorem, letting

	$\displaystyle\frac{1}{N}\sum_{n=1}^{N}A_{n}$	$\displaystyle=:a+\frac{\xi}{\sqrt{N}}\,,$
	$\displaystyle\frac{1}{N}\sum_{n=1}^{N}B_{n}$	$\displaystyle=:b+\frac{\eta}{\sqrt{N}}\,,$

where $a:=\mathbb{E}[A_{n}]$ , $b:=\mathbb{E}[B_{n}]$ , and $\xi\sim\mathcal{N}(0,\sigma_{a}^{2}),\ \eta\sim\mathcal{N}(0,\sigma_{b}^{2})$ follow the normal distribution with $\sigma_{a}^{2}:=\mathrm{Var}[A_{n}]$ , $\sigma_{b}^{2}:=\mathrm{Var}[B_{n}]$ , we have the expression (B.1) for $\bar{h}_{\nu,N}$ written as

(B.2)

\bar{h}_{\nu,N}=\frac{a+\xi/\sqrt{N}}{b+\eta/\sqrt{N}}\,.

To avoid the difficulty caused by the singularity with a probable zero value at the denominator in (B.2), we introduce the auxiliary function

(B.3)

f(x,y)=\frac{a+x}{\epsilon+g(b+y-\epsilon)}\,,

with a parameter $\epsilon\in(0,b/2)$ and a regularized piecewise function $g$ defined by

g(z)=\left\{\begin{aligned} &\quad 0\,,\ z<-\epsilon\,,\\ &\frac{(z+\epsilon)^{3}}{6\epsilon^{2}}\,,\ -\epsilon<z<0\,,\\ &\frac{(\epsilon-z)^{3}}{6\epsilon^{2}}+z\,,\ 0<z<\epsilon\,,\\ &\quad z\,,\ z>\epsilon\,.\\ \end{aligned}\right.

We have $g(z)\geq 0$ for all $z\in\mathbb{R}$ . Let $g_{\epsilon}(z):=\epsilon+g(z-\epsilon)$ , the function $f(x,y)$ defined in (B.3) can be written as

f(x,y)=\frac{a+x}{g_{\epsilon}(b+y)}\,.

A schematic plot of function $g_{\epsilon}(z)$ and a comparison between the functions $1/(b+y)$ and $1/g_{\epsilon}(b+y)$ is shown in Figure B.1.

If we have $y+b-\epsilon>0$ , then the two functions $1/(b+y)$ and $1/g_{\epsilon}(b+y)$ are identical. Given that the random variable $\eta\sim\mathcal{N}(0,\sigma_{b}^{2})$ and that $\epsilon<b/2$ , we have

		$\displaystyle\mathrm{Pr}\big{(}b+\frac{\eta}{\sqrt{N}}-\epsilon\leq 0\big{)}$
	$\displaystyle=$	$\displaystyle\mathrm{Pr}\big{(}\eta\leq-\sqrt{N}(b-\epsilon)\big{)}$
	$\displaystyle=$	$\displaystyle\big{\{}\eta\sim\mathcal{N}(0,\sigma_{b}^{2})\Leftrightarrow\eta=\sigma_{b}\,\omega,\textup{ where }\omega\sim\mathcal{N}(0,1)\big{\}}$
	$\displaystyle=$	$\displaystyle\mathrm{Pr}\big{(}\omega\leq-\frac{\sqrt{N}(b-\epsilon)}{\sigma_{b}}\big{)}$
	$\displaystyle=$	$\displaystyle\Phi\big{(}-\frac{\sqrt{N}(b-\epsilon)}{\sigma_{b}}\big{)}\,,$

where $\Phi$ is the cumulative distribution function of standard normal distribution. Recalling from Table 4.1 for our numerical experiments with inverse temperature $\beta=2$ , by employing a sample size $N=2^{28}$ to approximate the partition function $\bar{Z}(\beta)$ we obtain a 95% relative confidence interval $CI\approx 8\%$ , which gives

CI=1.96\frac{\sigma_{b}}{\sqrt{N}b}\approx 0.08\Longrightarrow\frac{\sqrt{N}b}{\sigma_{b}}\approx 24.5\,.

As the parameter $\epsilon$ satisfies $0<\epsilon<b/2$ , we obtain

\Phi(-24.5)\approx\Phi\big{(}-\frac{\sqrt{N}b}{\sigma_{b}}\big{)}<\Phi\big{(}-\frac{\sqrt{N}(b-\epsilon)}{\sigma_{b}}\big{)}<\Phi\big{(}-\frac{\sqrt{N}b}{2\,\sigma_{b}}\big{)}\approx\Phi(-12.3)\,,

which gives a tiny probability for a non-zero difference between $1/(b+\frac{\eta}{\sqrt{N}})$ and $1/g_{\epsilon}(b+\frac{\eta}{\sqrt{N}})$ , since the probability $\mathrm{Pr}\big{(}\eta\leq-\sqrt{N}(b-\epsilon)\big{)}$ decreases fast as $N$ increases.

Further we obtain

		$\displaystyle\frac{\partial f}{\partial x}=\frac{1}{\epsilon+g(b+y-\epsilon)},\quad\frac{\partial f}{\partial y}=-\frac{(a+x)(g^{\prime}(b+y-\epsilon))}{(\epsilon+g(b+y-\epsilon))^{2}}\,,$
		$\displaystyle\frac{\partial^{2}f}{\partial x^{2}}=0,\quad\frac{\partial^{2}f}{\partial x\partial y}=-\frac{g^{\prime}(b+y-\epsilon)}{(\epsilon+g(b+y-\epsilon))^{2}}\,,$
		$\displaystyle\frac{\partial^{2}f}{\partial y^{2}}=-\frac{(a+x)\big{(}g^{\prime\prime}(b+y-\epsilon)(\epsilon+g(b+y-\epsilon))+2\,g^{\prime}(b+y-\epsilon)^{2}\big{)}}{(\epsilon+g(b+y-\epsilon))^{3}}\,.$

and the third derivatives of $f$ are bounded by $\mathcal{O}(\frac{1}{\epsilon^{4}})$ . Applying Taylor’s theorem we have

	$\displaystyle f(x,y)=$	$\displaystyle\,f(0,0)+f^{\prime}_{x}(0,0)\,x+f^{\prime}_{y}(0,0)\,y+f^{\prime\prime}_{xx}(0,0)\,\frac{x^{2}}{2}+f^{\prime\prime}_{yy}(0,0)\,\frac{y^{2}}{2}+f^{\prime\prime}_{xy}(0,0)\,xy$
		$\displaystyle+\mathcal{O}\big{(}\frac{1}{\epsilon^{4}}(\|x\|^{3}+\|y\|^{3})\big{)}\,.$

With the condition $0<2\epsilon<b$ and substituting respectively $x$ and $y$ with $\xi/\sqrt{N}$ and $\eta/\sqrt{N}$ , based on (B.2) we further obtain,

(B.4)

\bar{h}_{\nu,N}\simeq\frac{a}{b}+\frac{1}{b}\frac{\xi}{\sqrt{N}}-\frac{a}{b^{2}}\frac{\eta}{\sqrt{N}}-\frac{\xi\eta}{b^{2}N}+\frac{a\eta^{2}}{b^{3}N}+\mathcal{O}\Big{(}\frac{1}{\epsilon^{4}}\big{(}\frac{|\xi|^{3}}{N^{\frac{3}{2}}}+\frac{|\eta|^{3}}{N^{\frac{3}{2}}}\big{)}\Big{)}.

Using the convexity of the functions $\xi^{3}$ and $\eta^{3}$ , we apply Jensen’s inequality to obtain

\Big{(}\mathbb{E}\Big{[}\frac{|\xi|^{3}}{N^{\frac{3}{2}}}\Big{]}\Big{)}^{\frac{4}{3}}\leq\mathbb{E}\Big{[}\frac{|\xi|^{4}}{N^{2}}\Big{]}=\mathcal{O}\big{(}\frac{1}{N^{2}}\big{)}\,,\quad\Big{(}\mathbb{E}\Big{[}\frac{|\eta|^{3}}{N^{\frac{3}{2}}}\Big{]}\Big{)}^{\frac{4}{3}}\leq\mathbb{E}\Big{[}\frac{|\eta|^{4}}{N^{2}}\Big{]}=\mathcal{O}\big{(}\frac{1}{N^{2}}\big{)}\,,

which yields

\mathbb{E}\Big{[}\frac{|\xi|^{3}+|\eta|^{3}}{N^{\frac{3}{2}}}\Big{]}=\mathcal{O}\big{(}\frac{1}{N^{\frac{3}{2}}}\big{)}\,.

Therefore, taking the expected value of $\bar{h}_{\nu,N}$ based on (B.4), we obtain

(B.5)

\mathbb{E}[\bar{h}_{\nu,N}]=\frac{a}{b}-\mathbb{E}\big{[}\frac{\xi\eta}{ab}\big{]}\frac{a}{b}\frac{1}{N}+\mathbb{E}\big{[}\frac{\eta^{2}}{b^{2}}\big{]}\frac{a}{b}\frac{1}{N}+\mathcal{O}\big{(}\frac{1}{\epsilon^{4}N^{\frac{3}{2}}}\big{)}\,.

Using specifically the covariance between $A_{n}$ and $B_{n}$ ,

	$\displaystyle\mathbb{E}[\xi^{2}]$	$\displaystyle=\mathbb{E}[(A_{n}-a)^{2}]=\sigma_{a}^{2}\,,$
	$\displaystyle\mathbb{E}[\eta^{2}]$	$\displaystyle=\mathbb{E}[(B_{n}-b)^{2}]=\sigma_{b}^{2}\,,$
	$\displaystyle\mathbb{E}[\xi\eta]$	$\displaystyle=\mathbb{E}[(A_{n}-a)(B_{n}-b)]=:\sigma_{ab}^{2}\,,$

we further obtain

(B.6)	$\displaystyle\mathbb{E}[(h_{\nu}-\bar{h}_{\nu,N})^{2}]$	$\displaystyle=\mathbb{E}\Big{[}\Big{(}\frac{\xi}{a\sqrt{N}}\frac{a}{b}-\frac{\eta}{b\sqrt{N}}\frac{a}{b}+\mathcal{O}\big{(}\frac{1}{\epsilon^{3}}(\frac{\xi^{2}}{N}+\frac{\eta^{2}}{N})\big{)}\Big{)}^{2}\Big{]}$
		$\displaystyle=\big{(}\frac{a}{b}\big{)}^{2}\frac{1}{N}\big{(}\frac{\sigma_{a}^{2}}{a^{2}}+\frac{\sigma_{b}^{2}}{b^{2}}-2\frac{\sigma_{ab}^{2}}{ab}\big{)}+\mathcal{O}\big{(}\frac{1}{\epsilon^{3}N^{\frac{3}{2}}}\big{)}$
		$\displaystyle=\big{\{}\,\textup{letting }\,\frac{\sigma_{a}^{2}}{a^{2}}+\frac{\sigma_{b}^{2}}{b^{2}}-\frac{2\sigma_{ab}^{2}}{ab}=:\sigma_{h}^{2}\big{\}}$
		$\displaystyle=\frac{1}{N}\big{(}\frac{a}{b}\big{)}^{2}\sigma_{h}^{2}+\mathcal{O}\big{(}\frac{1}{\epsilon^{3}N^{\frac{3}{2}}}\big{)}.$

According to the Berry-Esseen Theorem [16], for independent identically distributed random variables $X_{1}$ , $X_{2}$ , … with $\mathbb{E}[X_{1}]=0$ , $E[X_{1}^{2}]=\sigma^{2}$ , $\mathbb{E}[|X_{1}|^{3}]=\rho<+\infty$ , let $Y_{n}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ , then the cumulative distribution function (c.d.f.) $F_{n}$ of the normalized mean $\frac{\sqrt{n}Y_{n}}{\sigma}$ asymptotically approaches the c.d.f. of a standard normal distribution $\Phi$ as sample size $n$ increases. Specifically, $F_{n}$ satisfies

\sup_{n,x}|F_{n}(x)-\Phi(x)|\leq\frac{C\rho}{\sigma^{3}\sqrt{n}},\quad C<\frac{1}{2}\,.

Hence the stochastic variables $\xi$ and $\eta$ follow asymptotically normal distribution, with the difference between their c.d.f and the c.d.f. of their corresponding normal distributions bounded by $\mathcal{O}\big{(}\frac{1}{\sqrt{N}}\big{)}$ , where $N$ is the sample size for our mean-field estimator $\bar{h}_{\nu,N}$ . Practically, we also test with a series of sample sizes $N=2^{18}$ , $N=2^{20}$ , $N=2^{22}$ , and $N=2^{24}$ . For each fixed sample size value, we generate $M=256$ replicas of independent estimators $\{\tilde{Z}_{N}^{(m)}\}_{m=1}^{M}$ , each obtained by multiplying the partition function estimator ${\bar{Z}}_{N}^{(m)}$ with the constant factor $n!(2\pi\beta)^{\frac{3n}{2}}$ . By fitting a normal distribution of the obtained estimators for different $N$ values, we indeed observe the empirical standard deviation $\hat{\sigma}$ decreases with a factor around $1/2$ as the sample size $N$ increases with a factor $4$ , as can be observed in Figure B.2.

Moreover, the two subfigures corresponding to $N=2^{22}$ and $N=2^{24}$ show that out of the $M=256$ independent replicas, no observations with a negative value of $\tilde{Z}_{N}$ are present, validating our previous assertion that $\mathrm{Pr}\Big{(}\big{|}\frac{1}{\big{(}b+\frac{\eta}{\sqrt{N}}\big{)}}-\frac{1}{g_{\epsilon}\big{(}b+\frac{\eta}{\sqrt{N}}\big{)}}\big{|}>0\Big{)}$ is negligible with a sufficiently large sample size $N$ . Particularly in our numerical experiments, no bias related to the $\epsilon$ -regularization (B.3) is introduced.

In order to check that the conditions for applying the Berry-Esseen theorem are satisfied and to verify that based on our sample set of size $N$ , our estimator for the mean-field energy $\bar{h}_{\nu,N}\simeq\frac{a}{b}+\frac{1}{\sqrt{N}}\big{(}\frac{1}{b}\xi-\frac{a}{b^{2}}\eta\big{)}$ is approximately normally distributed, we also numerically evaluate the sample central moments of our estimators of $\mathcal{Z}_{\nu}$ and $\bar{h}_{\nu}$ . As can be observed from Figure B.3, the moments gradually stabilize as the sample size $M_{x}$ increases, and remain bounded up to the $4$ -th central moment. Combining the above reasoning with (B.5) and (B.6), we obtain

(B.7)

\bar{h}_{\nu,N}\simeq\frac{a}{b}\Big{(}1\pm\frac{\sigma_{h}}{\sqrt{N}}+\mathcal{O}\big{(}\frac{1}{N}\big{)}\Big{)}\,.

In practice, we use the sample mean and covariance values

		$\displaystyle\bar{a}=\frac{1}{N}\sum_{n=1}^{N}A_{n}\,,\quad\hat{\sigma}_{a}^{2}=\frac{1}{N-1}\sum_{n=1}^{N}(A_{n}-\bar{a})^{2}\,,$
		$\displaystyle\bar{b}=\frac{1}{N}\sum_{n=1}^{N}B_{n}\,,\quad\hat{\sigma}_{b}^{2}=\frac{1}{N-1}\sum_{n=1}^{N}(B_{n}-\bar{b})^{2}\,,$
		$\displaystyle\hat{\sigma}_{ab}^{2}=\frac{1}{N-1}\sum_{n=1}^{N}(A_{n}-\bar{a})(B_{n}-\bar{b})\,,$
		$\displaystyle\hat{\sigma}_{h}^{2}=\frac{\hat{\sigma}_{a}^{2}}{\bar{a}^{2}}+\frac{\hat{\sigma}_{b}^{2}}{\bar{b}^{2}}-\frac{2\hat{\sigma}_{ab}^{2}}{\bar{a}\bar{b}}\,,$

to obtain an empirical $95\%$ confidence interval for the mean-field estimator $\bar{h}_{\nu,N}$ by

\Big{(}\frac{\bar{a}}{\bar{b}}\big{(}1-1.96\frac{\hat{\sigma}_{h}}{\sqrt{N}}\big{)}\,,\ \frac{\bar{a}}{\bar{b}}\big{(}1+1.96\frac{\hat{\sigma}_{h}}{\sqrt{N}}\big{)}\Big{)}.

To ensure the reliability of the above statistical method (Method 1), we also implement an alternative statistical method (Method 2). For simplicity, the specific ideas for the Method 2 can be explained with the following example: if we want to employ a total sample size $M_{x}=2^{26}$ , we first obtain $M_{1}=2^{10}$ independent estimators of the mean-field energy $\{\bar{h}_{\nu}^{(i)}\}_{i=1}^{M_{1}}$ , and each of these estimators $\bar{h}_{\nu}^{(i)}$ is evaluated based on the independent sample sets $\{\mathbf{x}_{0}^{(i,j)},\mathbf{B}^{(i,j)}\}_{j=1}^{M_{2}}$ , where $M_{2}=2^{16}$ so that the total sample size $M_{x}=M_{1}\times M_{2}$ . By applying the Central Limit Theorem and fitting the obtained data $\{\bar{h}_{\nu}^{(i)}\}_{i=1}^{M_{1}}$ with Normal distribution, we obtain the statistical confidence intervals for $\bar{h}_{\nu}$ , and they are in good consistency with the confidence intervals obtained with Method 1 for sufficiently large $M_{1}$ and $M_{2}$ . Additional details for this implementation, along with sample code, are available on our open-access GitHub repository at https://github.com/XinHuang2022/Path_Integral_MD_mean_field.git.

Appendix C Particles with spin

When spin is included a transposition $T_{ij}$ requires the coordinates for the particle $i$ , in particular $(x_{i}(t),y_{i}(t),z_{i}(t),s_{i})$ , and particle $j$ to have the same spin, $s_{i}=s_{j}$ . If $s_{i}\neq s_{j}$ thev particles $i$ and $j$ are treated as distinguishable, since the Hamiltonian $\widehat{H}$ does not depend on the spin. We obtain two subsets of indistinguishable particles, one with spin $1/2$ and one with spin $-1/2$ , and determine the partition function for the system as follows. Let $\mathbf{x}=(\mathbf{x}_{+},\mathbf{x}_{-})$ , where $\mathbf{x}_{+}=(x_{1},y_{1},z_{1},s_{1},\ldots,x_{n_{+}},y_{n_{+}},z_{n_{+}},s_{n_{+}})$ is the coordinate for particles with spin $s_{1}=s_{2}=\ldots=s_{n_{+}}=1/2$ and $\mathbf{x}_{-}=(x_{n_{+}+1},y_{n_{+}+1},z_{n_{+}+1},s_{n_{+}+1},\ldots,x_{n},y_{n},z_{n},s_{n})$ is the coordinate for particles with spin $-1/2$ . An allowed coordinate permutation can then be written $\mathbf{x}^{\sigma}=(\mathbf{x}_{+}^{\sigma^{+}},\mathbf{x}_{-}^{\sigma^{-}})$ where $\mathbf{x}_{+}^{\sigma^{+}}$ is the $\sigma^{+}$ permutation among the particles with spin $1/2$ and $\mathbf{x}_{-}^{\sigma^{-}}$ is the $\sigma^{-}$ permutation among the particles with spin $-1/2$ . The antisymmetric wave functions have the representation

\phi(\mathbf{x})=\sum_{\sigma^{+}\in\mathbf{S}^{+}}\sum_{\sigma^{-}\in\mathbf{S}^{-}}\frac{{\rm sgn}(\sigma^{+})}{n_{+}!}\frac{{\rm sgn}(\sigma^{-})}{n_{-}!}\tilde{\phi}(\mathbf{x}^{\sigma})\,,

for $\tilde{\phi}\in L^{2}(\mathbb{R}^{3n})$ , which yields the partition function

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{3n}}\sum_{\sigma^{+}\in\mathbf{S}^{+}}\sum_{\sigma^{-}\in\mathbf{S}^{-}}\frac{{\rm sgn}(\sigma^{+})}{n_{+}!}\frac{{\rm sgn}(\sigma^{-})}{n_{-}!}\mathbb{E}[e^{-\int_{0}^{\beta}V(\mathbf{x}_{t},\mathbf{X}){\mathrm{d}}t}\delta(\mathbf{x}_{\beta}-\mathbf{x}_{0}^{(\sigma^{+},\sigma^{-})})]{\mathrm{d}}\mathbf{x}_{0}\,,

where $n_{\pm}$ particles have spin $\pm 1/2$ and the set of all coordinate permutations is denoted $\mathbf{S}^{\pm}$ . In the case of a separable potential

V(\mathbf{x},\mathbf{X})=\sum_{k=1}^{n_{+}}\tilde{V}(\mathbf{x}^{+}_{k},\mathbf{X})+\sum_{k=1+n_{+}}^{n}\tilde{V}(\mathbf{x}_{k}^{-},\mathbf{X})

the partition function becomes

\mathrm{Tr}\,(e^{-\beta H_{e}})=\int_{\mathbb{R}^{3n}}=\int_{\mathbb{R}^{3n}}\frac{{\rm det}\big{(}\mathcal{W}^{+}(\mathbf{x}^{+}_{0})\big{)}{\rm det}\big{(}\mathcal{W}^{-}(\mathbf{x}^{-}_{0})\big{)}}{n_{+}!n_{-}!(2\pi\beta)^{3n/2}}{\mathrm{d}}\mathbf{x}_{0}\,,

where $\mathcal{W}^{+}$ is the $n^{+}\times n^{+}$ matrix with components

\mathcal{W}_{k\ell}^{+}\big{(}\mathbf{x}^{+}\big{)}:=e^{-|\mathbf{x}_{k}^{+}(0)-\mathbf{x}_{\ell}^{+}(0)|^{2}/(2\beta)}e^{-\int_{0}^{\beta}\tilde{V}(\mathbf{B}_{k}(t)+(1-\frac{t}{\beta})\mathbf{x}_{k}^{+}(0)+\frac{t}{\beta}\mathbf{x}_{\ell}^{+}(0),\mathbf{X}){\mathrm{d}}t}\,,

and similarly for $\mathcal{W}_{-}$ .

(A.4)			$\displaystyle\frac{\partial}{\partial\beta}\,\gamma_{k,\ell}(\beta,\mathbf{x})=\frac{\partial}{\partial\beta}\Big{(}\int_{0}^{1}\big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\beta\,{\mathrm{d}}s\Big{)}$
			$\displaystyle\begin{aligned} =&\int_{0}^{1}\big{(}\frac{1}{2}\,\frac{\partial}{\partial\beta}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\partial}{\partial\beta}\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\beta\,{\mathrm{d}}s\\ &+\int_{0}^{1}\big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\big{)}\,{\mathrm{d}}s\\ \end{aligned}$
			$\displaystyle\begin{aligned} =&\int_{0}^{1}\Big{(}\frac{1}{2}\,\bar{\mathbf{x}}_{k}^{\ell}(s)\cdot\frac{1}{\sqrt{\beta}}\bar{\mathbf{B}}_{k}(s)\\ &-\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\big{(}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{)}}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}^{3}}\cdot\frac{1}{2\sqrt{\beta}}\big{(}\bar{\mathbf{B}}_{k}(s)-\bar{\mathbf{B}}_{j}(s)\big{)}\Big{)}\beta\,{\mathrm{d}}s\\ &+\int_{0}^{1}\Big{(}\frac{1}{2}\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}\bar{\mathbf{x}}_{k}^{\ell}(s)-\bar{\mathbf{x}}_{j}^{\nu_{j}}(s)\big{\|}}\Big{)}\,{\mathrm{d}}s\\ \end{aligned}$
			$\displaystyle\begin{aligned} =&\int_{0}^{\beta}\Big{(}\frac{1}{2}\,{\mathbf{x}}_{k}^{\ell}(t)\cdot\frac{1}{\beta}{\mathbf{B}}_{k}(t)\\ &-\frac{\lambda}{2}\sum_{j=1,j\neq k}^{n}\frac{\big{(}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{)}}{\big{\|}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{\|}^{3}}\cdot\frac{1}{2\beta}\big{(}{\mathbf{B}}_{k}(t)-{\mathbf{B}}_{j}(t)\big{)}\Big{)}\,{\mathrm{d}}t\\ &+\int_{0}^{\beta}\Big{(}\frac{1}{2}\big{\|}{\mathbf{x}}_{k}^{\ell}(t)\big{\|}^{2}+\sum_{j=1,j\neq k}^{n}\frac{\lambda}{2}\,\frac{1}{\big{\|}{\mathbf{x}}_{k}^{\ell}(t)-{\mathbf{x}}_{j}^{\nu_{j}}(t)\big{\|}}\Big{)}\frac{1}{\beta}\,{\mathrm{d}}t\,.\\ \end{aligned}$

Path integral molecular dynamics approximations of quantum canonical observables

Abstract.

1. Background to approximations of quantum canonical observables

2. Path integrals

2.1. Distinguishable particles

2.2. Indistinguishable particles

3. The fermion partition functions based on determinants

3.1. Motivation for the determinant approximation (3.8)

3.1.1. Marginals of 𝐱σj​(0)\mathbf{x}_{\sigma_{j}}(0)

3.1.2. The perturbation formulation (3.10)

Remark 3.1 (Bias in mean-field approximations).

4. Numerical experiments

4.1. Non-interacting Fermi gas in a confining harmonic trap

4.2. Electrons with Coulomb interaction under a harmonic trap

5. Conclusions

Acknowledgment

References

Appendix A Supplementary for Numerical Implementation

A.1. The derivative of the matrix element 𝒲k​ℓ​(β,𝐱)\mathcal{W}_{k\ell}(\beta,\mathbf{x}) with respect to β\beta

A.2. The time-rescaling technique on the Brownian bridge process 𝐁​(t)\mathbf{B}(t)

Appendix B Statistical confidence interval of the approximation of mean-field h¯ν\bar{h}_{\nu}

Appendix C Particles with spin

3.1.1. Marginals of $\mathbf{x}_{\sigma_{j}}(0)$

A.1. The derivative of the matrix element $\mathcal{W}_{k\ell}(\beta,\mathbf{x})$ with respect to $\beta$

A.2. The time-rescaling technique on the Brownian bridge process $\mathbf{B}(t)$

Appendix B Statistical confidence interval of the approximation of mean-field $\bar{h}_{\nu}$