Motivation

The two particle density matrix is of particular relevance for the understanding of quantum-mechanical many body systems. It contains more information than can be extracted from diagonal density matrix elements only and can be used to study eigenstates, their energies and the thermal excitations in the system. It is important for the characterization of the electronic excitations in the hydrogen plasma. In principle, one can determine the different electronic states and calculate the occupation probabilities. These are important for the radiative properties of the plasma and allow one to calculate absorption and emission spectra. The purpose of the following calculation is to determine the electronic excitations from a first principles calculation, which can then enter subsequent plasma models as an alternative to results from chemical models, in which the occupation numbers are estimated from the chemical equilibrium of approximately known eigenstates.

Furthermore, one would like to have a way to distinguish between free and bound states. This information can be used to estimate dynamic properties such as the conductivity. However, it has been proven to be a challenging task to give a rigorous and also practical criterion for the definition of free and bound states. In our approach, we followed the ideas of Girardeau (1990), who suggested the use of the two particle density matrix. Its eigenvalues are related to the occupation number and its eigenvector are called natural orbitals. They are of particular importance for configuration interaction methods in quantum chemistry because they lead to the fastest convergence (Szabo and Ostlund, 1996).

The two particle density matrix as defined in Eq. 5.2 requires path integral simulations with two open paths, which can be either of the same or of different particle type. One introduces relative and center-of-mass coordinates, ${\bf r}={\bf r}_1-{\bf r}_2$ and ${\bf r}_{{\rm cm}} = [m_1 {\bf r}_1 + m_2 {\bf r}_2 ] / (m_1+m_2)$. In a translationally invariant system, the pair matrix depends only on the separation of ${\bf r}_{{\rm cm}}-{\bf r}'_{{\rm cm}}$,

$$\rho^{[2]}({\bf r},{\bf r}',{\bf r}_{{\rm cm}},{\bf r}'_{{\r... ...\bf r},{\bf r}',{\bf r}_{{\rm cm}}-{\bf r}'_{{\rm cm}}) \quad,$$ (220)

and can be decomposed in its Fourier components (see Lebowitz et al. (1992) and Macris and Martin (1994)),

$$\rho^{[2]}_\mathbf{k}({\bf r},{\bf r}') = \int_{V^{\!\!\!\!\... ...m cm}}} \rho^{[2]}({\bf r},{\bf r}',{\bf r}_{{\rm cm}}) \quad.$$ (221)

If one averages over all momentum states one can write,

$$\rho^{[2]}({\bf r},{\bf r}') = \rho^{2}(r,r',\cos \theta) = ... ...2 l +1} \, \rho^{[2]}_l(r,r') \, {\rm P_l}(\cos \theta) \quad,$$ (222)

where $\theta$ is the angle between ${\bf r}$ and ${\bf r}'$. $\rho^{[2]}_l(r_,r')$ is a symmetric matrix with real and positive eigenvalues $n_{li}$. The eigenvectors $\phi_{li}(r)$ are the natural orbitals corresponding to a pair of particles in the angular momentum state $l$,

$$\rho^{[2]}_l(r_,r') = \sum_i \; n_{li} \; \phi_{li}(r) \; \phi^*_{li}(r') .$$ (223)

The eigenvalues $n_{li}$ indicate occupation probability of the $i$th orbital with angular momentum $l$ and energy $E_{li}$ given by,

$$n_{li} = \frac{1}{Z} \, e^{-\beta E_{li}} \quad,$$ (224)

$$Z = \sum_{li} e^{-\beta E_{li}} \quad.$$ (225)

It should be noted that the orbitals $\phi_{li}(r)$ and the energies $E_{li}$ correspond to many-particle states and are not the same as in single particle theory. However, the terminology used here is very similar.

In this approach, we did not make a distinction between the discrete spectrum of bound states and the continuum of free states. Both are included in the sum of Eq. 5.16. Since we applied this method to PIMC simulations in a finite volume, free states will have a discrete spectrum as well. In the simulation, one calculates the two particle density matrix from the following by average,

$$\rho^{[2]}_l(x,x') = \left< \: \delta(r-x) \: \delta(r'-x') \: {\rm P_l}(\cos \theta) \: \right> \quad.$$ (226)

In the MC method with open path, one can calculate the relative occupation numbers $n_{li}$. The Boltzmann factor, $e^{-\beta E_{li}}$, and therefore the orbital energies $E_{li}$ are not directly available because the partition function $Z$ cannot easily be calculated directly in any MC. However, one can, in principle, use the internal energy $E$ from a simulation with closed paths and sets,

$$E = \sum_{li} \; E_{li} \; n_{li} \quad,$$ (227)

which is correct in the mean field approximation. From Eq. 5.17, it follows,

$$\ln Z = -\beta E - \sum_{li} n_{li} \ln n_{li} \quad,$$ (228)

which is equivalent to the expression for the free energy $F=E-TS$. This means $Z$ can be calculated from the occupation numbers $n_{li}$ generated by an off-diagonal simulation and from the single particle energy $E$ taken from diagonal density matrix calculation. Then the energies of the orbital $E_{li}$ follow from Eq. 5.17. The disadvantage of this method is that it requires rather accurate estimates for all occupation numbers including highly excited states. As a consequence, we were not able to determine the absolute occupation numbers and energies even for a simulation of a single hydrogen atom in a periodic box, which will be discussed below.

This off-diagonal sampling method can be applied to different regimes in hydrogen. One can determine the electronic excitations in atoms by opening an electron and a proton path. A pair of open proton paths can be used to study the formation of hydrogen molecules. The eigenvectors then correspond to different rotational and vibrational excitations. In the following sections, we show how this method can be applied in PIMC simulations and show preliminary results for the electronic excitations in hydrogen.