Energetics and partition function of H molecular ion
Full quantum statistics of the H ion is simulated at low densities using the path integral Monte Carlo approach. For the first time, the molecular total energy, partition function, free energy, entropy and heat capacity are evaluated in temperatures relevant for planetary atmospheric physics. Temperature and density dependent dissociation recombination reaction balance of the molecule and its fragments above K is described, and also, the density dependence of thermal ionization above K is demonstrated. We introduce a new well-behaving analytical model for the molecular partition function of the H ion for the temperature range below dissociation and fit the parameters to the energetics from our simulations. The approach presented here can be regarded as an extension of the traditional ab initio quantum chemistry beyond the Born–Oppenheimer approximation to description of nonadiabatic phenomena, and even further, account of nuclear quantum dynamics.
Subject headings:molecular data – stars: atmospheres
The H molecular ion has been a subject of a number of theoretical and experimental studies since its first experimental detection (Thomson 1911). Because of its rapid formation through the exothermic reaction ( eV)
the H ion is expected in any active environment containing molecular hydrogen (Neale and Tennyson 1995), and thus, it is encountered, e.g. in the studies of the giant planets (Lystrup et al. 2008; Koskinen et al. 2009).
In planetary atmospheric physics, importance of the H ion lies in its capability to act as a cooling agent via infrared radiation (Neale et al. 1996; Harris et al. 2004; Koskinen et al. 2007). The atmospheric models taking into account this cooling are commonly based on the high temperature molecular partition function of the H ion (Neale and Tennyson 1995). Evaluation of the partition function faces, however, a few challenges of which the first one is finding a good approximation to the infinite summation over all rovibrational quantum states with accurate enough energies (Neale and Tennyson 1995). This has usually been worked out with calculations of a finite number of states from, e.g. a semi-empirical potential energy surface (Dinelli et al. 1995).
The next challenge comes with the changing geometry of the H ion at finite temperature. The rovibrational model needs to be extended for calculations of correct energetics for the emerging linear geometry of the weakly bound molecule (Neale et al. 1996).
Finally, at about the same temperature, where linear geometries start contributing, the molecule may also dissociate to its fragments, and in fact, the equilibrium
needs to be considered above about K, the balance depending strongly on both the temperature and the density.
This brings forth two questions, at the least. First, how relevant it is to consider the molecular energetics and related partition function at temperatures where the molecule has dissociated and appears in form of fragments of the equilibrium reaction, Eq. (1), only. Secondly, the balance of the equilibrium reaction may be strongly affected, not only by the density, but also by the environment including the neutralizing negative counterparts of the positive H. Thus, the thermal dissociation–recombination balance above about K gives rise to problems, which have not been taken into account in this context, yet.
In this study, using the path integral quantum Monte Carlo (PIMC) method we have carried out the first simulations of the full quantum statistics of the H ion, described by Eq. (1), at low densities and high temperatures ranging from K up to about K. PIMC is the method to meet the above challenges: we need not make any approximations or restrictions in the summing over states, geometries or quantum description of dynamics. The finite temperature is inherent in the PIMC approach and the Coulomb many-body treatment of the particle interactions is exact. The PIMC method is computationally expensive, but feasible for small enough systems. (Ceperley 1995; Pierce and Manousakis 1999; Kwon and Whaley 1999; Knoll and Marx 2000; Cuervo and Roy 2006; Kylänpää and Rantala 2009).
The conventional quantum chemical ab initio description of the H ion emerges as the zero Kelvin extrapolate from the PIMC simulations as we have shown earlier (Kylänpää and Rantala 2010). There, we evaluated the differences between three models for the description of the nuclear dynamics: Born–Oppenheimer approximation, nuclei in thermal motion and nuclei with both thermal and quantum dynamics. At low temperatures the necessity of the fully quantum mechanical approach for all five particles was established.
In the next section we present the essential details of the Feynman path integral quantum statistical approach, numerical simulation method and the model of the H ion. In section 3 we present and analyze the energetics, partition function and other thermodynamic functions of the system using analytical forms where pertinent. In the last section the conclusions are given.
2. Method and Models
According to the Feynman path integral formulation of the quantum statistical mechanics (Feynman 1998) the partition function of interacting distinguishable particles is given by the trace of the density matrix as
where the action is taken over the path in imaginary time , where and is called the Trotter number. The trace implies a closed path ().
For simulation we use the pair approximation in the action (Storer 1968; Ceperley 1995) for the Coulomb interaction of charges. This is exact in the limit , however, chemical accuracy is reached with sufficiently large , i.e. small enough . Sampling in the configuration space in ensemble is carried out using the Metropolis algorithm (Metropolis et al. 1953) with bisection moves and displacement moves (Chakravarty et al. 1998). The total energy is calculated using the virial estimator (Herman et al. 1982).
The error estimate in the PIMC scheme is commonly given in powers of the imaginary time time-step (Ceperley 1995). Therefore, in order to systematically determine the thermal effects on the system we have carried out all the simulations with , where denotes the unit of Hartree. Thus, the temperatures and the Trotter number are related by , where is the Boltzmann constant.
In the following we mainly use the atomic units, where the lengths, energies and masses are given in the units of Bohr radius (), Hartree () and free electron mass (), respectively. Thus, for the mass of the electrons we take and for the protons . Conversion of the units of energy is given by eV, and correspondingly, .
The statistical standard error of the mean (SEM) with SEM limits is used as an error estimate for the evaluated observables.
For the simulations we place one H ion, i.e. three protons and two electrons, into a cubic box and apply periodic boundary conditions and the minimum image principle. The simulations are performed in three different super cell (box) volumes: , and . These correspond to the mass densities of , and , respectively. The density has no essential effect at low , where dissociation rarely takes place. At higher , however, the finite density gives rise to the molecular recombination balancing the more frequent dissociation.
It should be pointed out that application of the minimum image principle with only one molecular ion in the periodic super cell may both give rise to the finite-size effects and also disregard higher density distribution effects, i.e. fragments of several ions in the simulation box. Thus, the lower the density the better we minimize the finite-size effects, which in this work are negligible, if not absent. The zero density limit cannot be reached due to the finite . To avoid all these ambiguities we define our targets as molecular energetics, molecular partition function and other related molecular quantities. Therefore, in the following, we also exclude the trivial contribution from the center-of-mass thermal dynamics and energy to the molecular quantities.
The nuclear quantum dynamics, which was shown to be essential at low , turns out to be negligible at higher temperatures. It is included, however, to be consistent with the low temperature results and our earlier study. For more details about the model and a discussion about the here neglected contribution from the exchange interaction see Kylänpää and Rantala (2010).
3. Results and discussion
3.1. Overview of molecular energetics
In Fig. 1 the total energy (canonical ensemble internal energy) of the H ion and its fragments is shown as a function of temperature. The molecular energy does not include the center-of-mass translational kinetic energy . The data from simulations are given as circles, squares and triangles corresponding to the three densities. The PIMC data is also given in Tables 1 and 2.
The solid lines at K are fitted to analytical model forms but at higher temperatures lines are for guiding the eye, only. Our low temperature fit and analytical model, Eq. (7), is given as a blue dashed line and is discussed in the next sections in more detail. For comparison the energies from the fitted partition function of Neale and Tennyson (1995) is shown as black dots. These two do not manifest dissociation, and therefore, are not relevant at "higher ".
The horizontal dash-dotted lines show the zero Kelvin energies for the ion and its fragments in Eq.(1). One of these lines presents the energy for the "barrier to linearity", i.e. the energy needed for the transformation to linear molecular geometry. It is also seen to be roughly the barrier to dissociation within the considered molecular densities.
Above K the density dependence is clearly seen as varying composition of fragments. In the range from to K the changing dissociation recombination balance leads to distinctly different energetics, and above that, at our highest simulation temperatures the thermal ionization of hydrogen atoms starts contributing to the energy. However, it is worth pointing out that the temperature limits of these three ranges, i.e. K, K and above K, are subject to changes with larger variation of densities.
Above K in our lowest density case the thermal ionization of H atoms is evident, see Fig. 1, but for our higher density cases some K is needed to show first signs of ionization. Similar trend for the ionization is stated in Koskinen et al. (2010), although there the density is notably less than our lowest one.
Let us now consider the dissociation recombination reaction chain, Eq. (1), and the contributing fragments to the quantum statistical equilibrium trying to give an intuitive classical-like picture of the composition. With finite , instead of zero, we have finite , instead of infinite, that brings classical nature to the system the more, the higher the temperature. In other words, the partial decoherence in our five particle quantum system increases with increasing temperature, that enables us to distinguish the fragments as separate molecules and atoms in thermal equilibrium. Based on this interpretation, we show the total energy distribution in Fig. 2 from sampling the imaginary time paths at about K with in our highest simulation density.
We see three main peaks and by inspection of the energy distribution the first and the second can clearly be assigned to the rovibrationally excited H and HH, respectively. As expected, there are no rovibrational excitations available for HH, leading to the peak average position very close to . The fourth fragment, HH, can be identified as the small high-energy side shoulder of HH peak. With the interpretation of the area under the peak as the abundance of the fragment in the equilibrium we find this contribution to be much smaller than that of the others.
It is important to note, however, that the above illustration is dependent on the block averaging procedure, see the caption of Fig. 2. Pinning the energy data of each and every sample, i.e. choosing block of size one sample, would broaden the peaks in Fig. 2. At the opposite limit, all samples in one block, the single mean energy or ensemble average is the quantum statistical expectation value , see Table 2, Figs. 1 and 2, where the statistical uncertainty decreases with increasing simulation length.
3.2. Molecular partition function
To compare with the other published approaches for the molecular partition function based on single molecule quantum chemistry we start from the lowest temperature range from to K, where the molecule does not essentially dissociate, yet.
We present a low temperature H molecular partition function as a first approximation for the modeling of low density H ion containing atmospheres. Our aim is to find a simple analytical form, which can be accurately fitted to the energies from our simulations.
The partition function in terms of the Helmholtz free energy is written as
where , and the energy expectation value is straightforwardly derived from the partition function as
After solving the free energy from Eq. (3) as
we write and the energy expectation value may be written as
We find that a well-behaving function fitting perfectly into our simulation data,
allows analytical integration of Eq. (6) for or ,
Using the boundary condition for the molecular partition function with a nondegenerate ground state, or , we get in our model. However, inclusion of the contributions from the ground state spin degeneracy factor and the zero-point rotations would give , which would lead to , and thus, shift the function by a constant, only.
In the fit, in addition to the weights, we force the first derivative of the energy with respect to the temperature to be monotonically increasing up to K. The fit extrapolates the K energy to about above that of the para-H, i.e. it gives an excellent match within the statistical error estimate.
In Fig. 3 the function from Eq. (8) is shown in the range K — the behaviour of the model at higher is illustrated by the dashed line. Above K the three curves for different densities are obtained from those shown in Fig. 1 by numerical integration of Eq. (6) as
Neale and Tennyson (1995) have presented the partition function based on a semi-empirical potential energy surface, see Fig. 3. The overall shape is similar to the one of ours. However, the energy evaluated from their fit tends to be systematically lower than ours, although roughly within our 2 SEM error limits. Thus, the deviations are not visible in Fig. 1. For the partition function the main difference is in their zero reference, which also leads to . This difference in zero reference goes back to that of energetics: our in Eq. (7) at fits perfectly the zero Kelvin energy of para-H, as mentioned earlier. The zero Kelvin energy of the Neale and Tennyson (1995) behind their partition function is that of the spin forbidden state, i.e. . This in fact, leads to the divergence of at .
Our low temperature partition function Eq. (8) is close to complete. With the PIMC approach we implicitly include all of the quantum states in the system with correct weight without any approximations. This partition function is the best one for the modeling of the low density H ion containing atmospheres, at the moment. However, it is valid up to about K, only. As soon as the density dependence starts playing larger role more complex models are needed. Such models can be fitted to our PIMC data given in Tables 1 and 2.
3.3. Other thermodynamic functions
In Fig. 4 we show the Helmholtz free energy from Eqs. (5) and (8), combined. As expected, lower density or larger volume per molecule lowers the free energy due to the increasing entropic factor. Dissociation and the consequent fragments help in filling both the space and phase space more uniformly or in less localized manner.
This kind of decreasing order is seen more clearly in the increasing entropy, shown in Fig. 5. The entropy has been evaluated from
where the internal energy is . As expected, both the total energy (internal energy) and entropy reveal the dissociation similarly.
Finally, in Fig. 6 we present the heat capacity
where is taken from Eq. (7), which is valid at low temperatures, only.
Considering the goodness of our functional form for , it is very convincing to see the plateau at about corresponding to "saturation" of the contribution from the three rotational degrees of freedom. Thus, above K the rotational degrees of freedom obey the classical equipartition principle of energy. It is the last term in the functional form of Eq. (7), that gives the flexibility for such detailed description of the energetics.
It should be emphasized that the plateau is not artificially constructed to appear at , except for a restriction given for the first derivative of the total energy to be increasing. Thus, the analytical model we present, Eq. (7), is found to be exceptionally successful at low temperatures.
We have evaluated the temperature dependent quantum statistics of the five particle molecular ion H far beyond its dissociation temperature at about K. This is done with the path integral Monte Carlo (PIMC) method, which is basis set and trial wave function free approach and includes the Coulomb interactions exactly. Thus, we are able to extend the traditional ab initio quantum chemistry with full account of correlations to finite temperatures without any approximations, also including the nuclear thermal and quantum dynamics.
It is fair to admit, however, that PIMC is computationally heavy for good statistical accuracy and approximations are needed to solve the "Fermion sign problem" in cases where exchange interaction becomes essential.
The temperature dependent mixed state description of the H ion, the density dependent equilibrium dissociation recombination balance and the energetics has been evaluated for the first time. With the rising temperature the rovibrational excitations contribute to the energetics, as expected, whereas the electronic part remains in its ground state in the spirit of the Born–Oppenheimer approximation. At about K the fragments of the molecule, HH, HH and HH, start contributing. Therefore, H ion becomes less dominant, and eventually negligible in high enough .
We have shown how the partial decoherence in the mixed state can be used for interpretation of the fragment composition of the equilibrium reaction. Furthermore, we have evaluated explicitly the related partition function, free energy, entropy and heat capacity, all as functions of temperature. An accurate analytical functional forms are given for the temperatures below dissociation. We consider all these as major improvements to the earlier published studies, where dissociation has not been considered.
For financial support we thank the Academy of Finland, and for computational resources the facilities of Finnish IT Center for Science (CSC) and Material Sciences National Grid Infrastructure (M-grid, akaatti).
- Ceperley (1995) Ceperley, D. M. 1995, Rev. Mod. Phys, 67, 279
- Chakravarty et al. (1998) Chakravarty, C., Gordillo, M. C., and Ceperley, D. M. 1998, J. Chem. Phys., 109, 2123
- Cuervo and Roy (2006) Cuervo, J. E. and Roy, P.-N. 2006, J. Chem. Phys., 125, 124314
- Dinelli et al. (1995) Dinelli, B. M., Polyansky, O. L., and Tennyson, J. 1995, J. Chem. Phys., 103, 10433
- Feynman (1998) Feynman, R. P. 1998, Statistical Mechanics, Perseus Books, Reading, Massachusetts
- Harris et al. (2004) Harris, G. J., Lynas-Gray, A. E., Miller, S., and Tennyson, J. 2004, ApJ, 600, 1025
- Herman et al. (1982) Herman, M. F., Bruskin, E. J., and Berne, B. J. 1982, J. Chem. Phys., 76, 5150
- Knoll and Marx (2000) Knoll, L. and Marx, D. 2000, Europ. Phys. J. D, 10, 353
- Koskinen et al. (2009) Koskinen, T. T., Aylward, A. D., and Miller, S. 2009, Astrophys. J., 693, 868
- Koskinen et al. (2007) Koskinen, T. T., Aylward, A. D., Smith, C. G. A., and Miller, S. 2007, ApJ, 661, 515
- Koskinen et al. (2010) Koskinen, T. T., Yelle, R. V., Lavvas, P., and Lewis, N. K. 2010, ApJ, 723, 116
- Kutzelnigg and Jaquet (2006) Kutzelnigg, W. and Jaquet, R. 2006, Phil. Trans. R. Soc. A, 364, 2855
- Kwon and Whaley (1999) Kwon, Y. and Whaley, K. B. 1999, Phys. Rev. Lett., 83, 4108
- Kylänpää and Rantala (2009) Kylänpää, I. and Rantala, T. T. 2009, Phys. Rev. A, 80, 024504
- Kylänpää and Rantala (2010) Kylänpää, I. and Rantala, T. T. 2010, J. Chem. Phys., 133, 044312
- Lystrup et al. (2008) Lystrup, M. B., Miller, S., Russo, N. D., R. J. Vervack, J., and Stallard, T. 2008, ApJ, 677, 790
- Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E. 1953, J. Chem. Phys., 21, 1087
- Neale et al. (1996) Neale, L., Miller, S., and Tennyson, J. 1996, ApJ, 464, 516
- Neale and Tennyson (1995) Neale, L. and Tennyson, J. 1995, ApJ, 454, L169
- Pavanello and Adamowicz (2009) Pavanello, M. and Adamowicz, L. 2009, J. Chem. Phys., 130, 034104
- Pierce and Manousakis (1999) Pierce, M. and Manousakis, E. 1999, Phys. Rev. B, 59, 3802
- Storer (1968) Storer, R. G. 1968, J. Math. Phys., 9, 964
- Thomson (1911) Thomson, J. J. 1911, Philos. Mag., 21, 225