Accuracy of the Pair Density Matrix

Figure: Internal energy (left graph) and $2K+V$ (right graph) for an isolated hydrogen atom as a function of the number of time slices at constant temperature of $10\,000\,\rm K$.

The pair density matrices are calculated using the matrix squaring method described in section 2.3.1. They are stored in tables using the expansion formula Eq. 2.38 and then entered into the PIMC simulation program. The accuracy of these tables is crucial for all following results. Using the precomputed pair density matrices allows one to employ a much larger time step because one starts with a solution of the two-particle problem. Fig.  shows how accurate this method is. The internal energy of an isolated hydrogen atom at sufficiently temperature ( $T=10\,000\,\rm K$) in a large box ($L=26$) is compared with the exact groundstate energy of $-13.6\,\rm eV$. The temperature was chosen low enough so that excited states can be neglected e.g. the contribution to the energy from the occupation of first excited state is $7 \cdot 10^{-5}\,\rm eV$ at this temperature. Furthermore, it is tested whether the kinetic energy $K$ and the potential energy $V$ satisfy the virial theorem $2K+V=0$. If only diagonal action terms are considered in Eq. 2.38 one finds a rather slow convergence as function of the number of time slices (Fig. ). Eventually, the error goes to zero in the Trotter limit, Eq. 2.24, of an infinite number of slices. Using off-diagonal terms in the expansion formula, Eq. 2.38, improves the convergence significantly as shown in Fig. . One can use different orders to calculate the action $n_{\rm A}$ and the energy $n_{\rm E}$. The resulting accuracy from different orders is shown in Tabs.  and . It reveals that using order 2 or higher instead of order 1 in the action decreases the errors by almost one order of magnitude. This is an important observation because most many-particle simulations reported in this work had been performed with $n_{\rm A}=1$ and $n_{\rm E}=2$ (which was found to be sufficient for simulations of Helium particles) before this analysis was done. Because of the resulting inaccuracies, the estimated energies and pressures are slightly too high.

Table: Accuracy study of PIMC simulations of the isolated hydrogen atom using different orders in the expansion formula 2.38 for the action and the energy. The calculated $2K+V$ (exact value equals zero) and the deviation of the potential energy $V$ from the exact value of -27.2$\,$eV are listed from simulations at $T=10\,000\,\rm K$ using 100 time slices and $\tau^{-1}=10^{6}\rm K$.

	$2K+V\,\rm (eV)$				$V-V_{\rm exact}\,\rm (eV)$
energy	action order				action order
order	0	1	2	3	0	1	2	3
0	0.8(1)				1.74(3)
1		0.47(2)	-0.18(3)	-0.071(25)		-0.176(08)	0.074(11)	0.043(07)
2		0.56(2)	-0.07(3)	0.013(48)		-0.166(05)	0.063(09)	0.034(13)
3		0.61(4)	-0.06(3)	0.039(08)		-0.183(11)	0.056(11)	0.031(03)

Table: Accuracy study of PIMC simulations using different order in the action and energy expansion as in Tab.  but with a larger time step of $\tau^{-1}=0.5 \cdot 10^{6}\rm K$ and 50 time slices.

	$2K+V\,\rm (eV)$				$V-V_{\rm exact}\,\rm (eV)$
energy	action order				action order
order	0	1	2	3	0	1	2	3
0	1.24(7)				1.991(28)
1		0.46(4)	-0.12(5)			-0.074(12)	0.033(12)
2		0.51(4)	-0.07(2)			-0.067(13)	0.044(07)
3				0.03(2)				0.036(7)

Table: Accuracy analysis for the isolated hydrogen molecule for different time steps and temperatures using $n_{\rm A}=3$ and $n_{\rm E}=3$. In this analysis, the nuclei are classical and fixed at the bond length of $R=1.4008$. $2K+V$ (exact value equals zero) and the deviations from the exact binding potential energy per atom $-31.946\,\rm{eV}$ are listed.

	$2K+V\,\rm (eV)$			$V-V_{\rm exact}\,\rm (eV)$
	$\tau^{-1} \; (10^6 {\rm K})$			$\tau^{-1} \; (10^6 {\rm K})$
$T ({\rm K})$	0.5	1	2	0.5	1	2
7812	0.41(3)	0.16(3)	-0.06(8)	-0.422(10)	-0.112(8)	0.030(22)
3906	0.41(3)	0.19(3)	0.14(4)	-0.419(07)	-0.110(5)	-0.008(11)
1953	0.43(3)	0.16(2)	0.08(3)	-0.421(12)	-0.115(9)	-0.009(09)

Figure: Internal energy (right graph) and $2K+V$ (left graph) per atom for an isolated hydrogen molecule as a function of the nuclear separation $R$. The PIMC results were calculated at T=3906$\,$K using a time step of $\tau^{-1}=10^6\,\rm K$ and the orders in the action and energy expansion ($n_{\rm A}$,$n_{\rm E})$ as indicated. The dashed line in the left graph indicates zero of energy and in the right graph, the energy of one isolated atom.

These inaccuracies have also consequences for the isolated hydrogen molecule, which are important to study in order to determine the corrections of results from many particle simulations. The virial theorem for an isolated H$_2$ molecular with fixed nuclei at separation $R$ reads (Kolos and Wolniewicz, 1964; Steiner, 1976),

$$2K+V = - R \; \frac{dE}{dR} \quad.$$ (198)

This means one can calculate the kinetic and potential energy from $E(R)$ and its derivative using $E=K+V$. In Fig. , the exact results by Kolos and Wolniewicz (1964) for $E$ and $2K+V$ are compared with the PIMC calculations. It shows that for PIMC using $\tau^{-1}=10^6\,\rm K$, $n_{\rm A}=1$, and $n_{\rm E}=2$, the energy $E$ is too high by $0.25 \pm 0.05\rm\,eV$ and $2K+V$ is too large by $0.7\pm 0.1\,\rm eV$ per atom. The latter correction is particularly important because in PIMC simulations of Coulomb systems, we use the virial theorem, Eq. 3.55, to estimate the pressure. The correction to the pressure is equivalent to subtracting the pressure of an ideal H$_2$ gas at $T=5400\pm800\,\rm K$.

The accuracy of PIMC simulations of an isolated molecule is affected by the order in action and energy expansion as well as by the time step because it is a four-particle problem. Tab.  shows results calculated with $n_{\rm A}=n_{\rm E}=3$ orders. First, we studied the different temperatures and found no dependence on $T$, which means that contributions from electronic excited states are negligible. Furthermore, the comparison of different time steps shows a significant dependence. Using a time step $\tau^{-1}=2\cdot10^6\,\rm K$ allows one to calculate the energy with an accuracy of approximately $0.05\,$eV and $2K+V$ with an error of about $0.1\,$eV per atom. Using a smaller time step would bring the results in Fig.  closer to the exact results.

Most many-particle PIMC simulations discussed in the following sections have been performed using $n_{\rm A}=1$ and $n_{\rm E}=2$ in the action and energy expansion of the pair density matrices. The correction resulting from higher order terms will be estimated based on the following argument. Higher order off-diagonal terms are large for small separations of the two particles. Therefore, we expected the dominant corrections to come from pairs of protons and electron when both particles are close together. Therefore, we suggest to use the integral of the proton-electron pair correlation function up to a cut-off radius,

$$I_{\rm pe} = 4 \pi n \int_0^{r_c} \!\!\!\!dr \; r^2 \, g_{\rm pe} (r) \quad,$$ (199)

in order to estimate the corrections to pressure and energy,

$$\Delta E = c_1 \; I_{\rm pe}\quad,$$ (200)

$$\Delta p = c_2 \; I_{\rm pe} \; \frac{n}{3}\quad.$$ (201)

We estimated $r_c=1.4$ from studying the magnitude of the higher order terms and determined the coefficients $c_1$ and $c_2$ from the corrections for the isolated molecule discussed above: $c_1=-0.24 \pm 0.05\rm\,eV$ and $c_2=-0.7\pm 0.1\,\rm eV$. In the following discussion of the thermodynamic properties, it will be explicitly stated where this correction has been applied. It turns out that the corrections to the pressure are particularly important at low densities and temperatures, where the pressure becomes of the same order of magnitude as the correction. It should also be noted that the corrections hardly change the hugoniot curve discussed in section . The effect is significantly smaller than the differences in energy and pressure, which are relevant in this context.