^{†}Corresponding author. Email: hongzhang@scu.edu.cn
^{*}Project supported by the National Natural Science Foundation of China (Grant Nos. 11474207 and 11374217).
The path integral Monte Carlo (PIMC) method is employed to study the thermal properties of C_{70} with one, two, and three H_{2} molecules confined in the cage, respectively. The interaction energies and vibrationally averaged spatial distributions under different temperatures are calculated to evaluate the stabilities of (H_{2})_{ n}@C_{70} ( n = 1, 2, 3). The results show that (H_{2})_{2}@C_{70} is more stable than H_{2}@C_{70}. The interaction energy slowly changes in a large temperature range, so temperature has little effect on the stability of the system. For H_{2}@C_{70} and (H_{2})_{2}@C_{70}, the interaction energies keep negative; however, when three H_{2} molecules are in the cage, the interaction energy rapidly increases to a positive value. This implies that at most two H_{2} molecules can be trapped by C_{70}. With an increase of temperature, the peak of the spatial distribution gradually shifts away from the center of the cage, but the maximum distance from the center of H_{2} molecule to the cage center is much smaller than the average radius of C_{70}.
Endohedral fullerene complex, with small atoms or molecules confined inside an empty fullerene cage, has attracted increasing attention. Because of the interaction between the host fullerene cage and the guest atoms, endohedral fullerene complex not only has the properties of the host cage but it also has the properties of the guests, ^{[1, 2]} which makes it widely used in many fields, such as electronics, ^{[3– 5]} medicine, ^{[5– 7]} environment protection, ^{[8]} and renewable energy.^{[9]}
Since La@C_{82} was macroscopically prepared, ^{[10]} significant progress has been made in researching endohedral metallofullerenes. It has been found that a fullerene complex with noble gas atoms can be prepared at high temperature (650 ° C) and high pressure (3000 atm, 1 atm = 1.01325 × 10^{5} Pa).^{[11, 12]} Much theoretical work has also been done to study the stability of this kind of complex.^{[13– 16]} In recent years, a fullerene complex with one or multiple hydrogen molecules has become a popular topic in both experimental and theoretical fields. The amount of work in which various theoretical methods are used has been done with C_{60} to study the maximum number of H_{2} molecules that can be trapped by the cage. Turker and Erkoc^{[17]} reported that 24 hydrogen molecules could be inserted into C_{60} from semiempirical Austin Model 1 (AM1) calculations. Moreover, Ren et al.^{[18]} reported that C_{60} could encapsulate up to 25 hydrogen molecules on the basis of parametric method 3 (PM3) combined with densityfunctional theory (DFT). However, the methods they used were considered to be defective and their results were thought to be unreliable.^{[19, 20]} Tatiana Korona et al.^{[21]} used the density functional theory combined with the symmetry adapted perturbation theory (DFTSAPT) and found that only one H_{2}molecule could be stabilized inside C_{60}. This result is consistent with the fact that H_{2}@C_{60} has been prepared but not (H_{2})_{2}@C_{60}.
Recent syntheses of endohedral complexes of C_{70} with one or two hydrogen molecules have sparked the interest of researchers who wish to study the stability of one or multiple H_{2} molecules encapsulated into C_{70} and to predict the maximum number of H_{2} that can be trapped by C_{70}. Some quantum mechanical methods, such as DFT, secondorder Mø ller– Plesset theory (MP2), spincomponent scaled MP2 (SCSMP2), diffusion Monte Carlo (DMC) have been used to evaluate the stabilities of H_{2}@C_{70} and (H_{2})_{2}@C_{70} and they found that both H_{2}@C_{70} and (H_{2})_{2}@C_{70} are energetically stable.^{[22– 24]} Sebastianelli et al.^{[24]} calculated the vibration properties of one or two H_{2} molecules in a C_{70} cage. However, their evaluation of the stability and the study of the thermal properties were conducted under the condition of temperature T = 0. Therefore their study cannot reflect the effect of temperature on the properties of the system. Since temperature is a parameter of the density matrix, the path integral Monte Carlo (PIMC) method is an accurate tool that can be used to compute quantum manybody properties at finite temperatures.^{[25– 27]} In this paper we employ the PIMC method to calculate the thermal properties of (H_{2})_{n}@C_{70} (n = 1, 2, 3).
The rest of this paper is organized as follows. In Section 2, we briefly introduce part of the PIMC theory relating to our study and the details of our simulation. In Section 3, we give our results and discussion, analyzing the temperaturedependent interaction energies of (H_{2})_{n}@C_{70} (n = 1, 2, 3) obtained from PIMC calculations and we will study the movement of H_{2} molecules in the C_{70} cage through onedimensional (1D) and threedimensional (3D) probability distribution functions (PDFs). Our conclusions are given in Section 4.
In the PIMC method, all observations can be expressed as statistics relating to the sampled paths. In the coordinate image, when a system is in a state of thermal equilibrium, the value of an operator O can be written as
where the partition function is z = ∫ drρ (R, R; β ), with ρ (R, R′ ; β ) denoting the density matrix.
In quantum mechanics, the density matrix is a basic physical quantity. Theoretically, all of the equilibrium properties of a quantum system can be obtained from the density matrix, which can be written as^{[25]}
where β is the inverse temperature and β = 1/K_{B}T, with K_{B} being the Boltzmann constant and T being the temperature, and H is the Hamiltonian. Since the density matrix of the system is not just about the Hamiltonian but also about the system temperature, the PIMC method can calculate the properties of systems under different temperatures.
The above expression can be expanded into a path of M steps in imaginary time, each step of duration τ = β /M is expressed as^{[25, 27]}
where {R, R_{1}, R_{2}, … , R_{M− 1}, R′ } denotes the path in position space. Hence the density matrix at temperature T can be expressed as the density matrix at higher temperature MT. When M is a limited integer, the path is discrete; however, when M → ∞ , the path is continuous.
Because multiple H_{2} molecules are involved in our system, the exchange between these identical particles must be taken into account. Since H_{2} molecules are bosons, we only consider the exchange occurring in a Bose system. For identical particles, the density matrix is^{[24]}
where P is a permutation, N is the number of the particles, and (P_{1}, P_{2}, … , P_{N}) is an allpermutations of (1, 2, … , N). The system with N identical particles has N! different permutations.
PDFs are used to describe the vibrationally averaged spatial distribution of H_{2} molecules, ^{[28]} which are defined as
Equation (5) is a 1D probability function and equation (6) is a 3D probability function, which can be calculated by the PIMC method, (x, y, z) are the Cartesian coordinates of the H_{2} molecule, 〈 · · · 〉 denotes the configurational average over the sampled paths, and R is the distance from the cage center to the center of caged H_{2} molecule. The PDFs are different under different temperatures and, when temperature is considered, the PDFs can be written as P(R; T) and P(x, y, z; T).
Because of the encapsulation a small number of H_{2} molecules can cause negligible distortions to the fullerene geometries, ^{[29]} we treat fullerene C_{70} as being rigid. In the case of hydrogen molecules, we treat them as spatial particles, which is proven to be accurate by calculating the properties of confined H_{2} in isotropic environments.^{[30, 31]} In our study, the C_{70} cage is considered to be fixed, our study only involves the movement of H_{2} molecules. The (H_{2})_{n}– C_{70} (n = 1, 2, 3) interaction is described by the sum of the pair potentials between H_{2} molecule and each of the 70 carbon atoms, as well as the intermolecular pair potentials between each H_{2} molecules when n > 1. For H_{2}– C pair potential, we use 6– 12 LennardJones (LJ) potential (6– 12 LJ potential).^{[32]} In the case of H_{2}– H_{2} pair potential, Silver– Goldman potential (SG potential)^{[33]} is employed. In our program, we use linear grids and the number of grids is 250. For carbon atoms, which are considered to be classical particles, the quantumness parameter λ is taken as 0 K· Å ^{2}. For H_{2} molecules, λ is taken to be 12.127 K· Å ^{2}. The number of partial waves is set to be 100. nsquare, which is the total number of squarings to reach the lowest temperature and is set to be 14. We use the bisection method to sample the paths and ensure that the acceptance for the moves reaches about 50%.
In this paper, temperaturedependent interaction energies are used to evaluate the stabilization of the endohedral fullerene complexes (H_{2})_{n}@C_{70} (n = 1, 2, 3), which are shown in Fig. 1. Table 1 shows our results at 1 K, which are compared with the results from three other quantummechanical calculations under the condition of temperature T = 0.
At 1 K, the interaction energies of H_{2}@C_{70}, (H_{2})_{2}@C_{70}, and (H_{2})_{3}@C_{70} calculated by PIMC method are − 1234.53 cm^{− 1}, − 1648.07 cm^{− 1}, and 2670.75 cm^{− 1}, respectively. This result is closest to the values obtained by DMC calculation (T = 0). The maximum difference between them is about 11%. The reason for this difference is that the interaction potentials used in DMC calculation and in this work are different. In this work, we use the SG potential to describe the interaction of H_{2}– H_{2}. In the DMC method, the H_{2}– H_{2} potential is described by the ab initio fourdimensional (4D) intermolecular potential energy surface (PES). For H_{2}– C pair potential, we use 6– 12 LJ potential, which is a standard twosite H_{2}– C pair potential; however, the DMC method uses the threesite H_{2}– C pair potential, which can be used to study the translation– rotation excitations of the system.
Since the interaction energies for one and two H_{2} molecules caged in C_{70} are negative, it can be inferred that both H_{2}@C_{70} and (H_{2})_{2}@C_{70} are energetically stable. Furthermore, the interaction energy of (H_{2})_{2}@C_{70} is greater than that of H_{2}@C_{70}, about a 34% increase (in absolute value), implying that the former is more stable than the latter. However, when three H_{2} molecules are encapsulated in C_{70}, the interaction energy rapidly jump to a high positive value (2670.75 cm^{− 1}). The reason why C_{70} cannot trap three H_{2} molecules may be that the C_{70} room is too crowded for three H_{2} molecules to move and each H_{2} molecule will be repelled by the other two H_{2} molecules.
From Fig. 1 we can see that for each of the complexes (H_{2})_{n}@C_{70} (n = 1, 2), the interaction energy increases gradually with the increase of temperature, which happens because of a stronger vibration of H_{2} molecules at a higher temperature. Meanwhile, in a wide temperature range (from 1 K to 300 K), the interaction energy remains negative for H_{2}@C_{70} and (H_{2})_{2}@C_{70}, so temperature has little effect on interaction energy. The result calculated by the PIMC method shows that C_{70} can trap two H_{2} molecules at most and is in agreement with the syntheses of H_{2}@C_{70} and (H_{2})_{2}@C_{70} but not (H_{2})_{3}@C_{70}.
The 1D and 3D PDFs are calculated to analyze the spatial distributions of one and two H_{2} molecules in the C_{70} cage. Figure 2 displays the 3D PDF, P(x, y, z), under 1 K from the PIMC calculations for (H_{2})_{n}@C_{70} (n = 1, 2). This figure can vividly show the spatial distribution of H_{2} molecules in C_{70}. For n = 1, the P(x, y, z) is eggshaped and occupies the center of the cage, meanwhile it extends along the C_{5} axis of C_{70}. For n = 2, the P(x, y, z) consisting of two distinct lobes, which are symmetrically distributed on both sides of the cage center, that the C_{5} axis of C_{70} passes through. What is noteworthy is that the bottoms of the lobes are flat, proving that the two H_{2} molecules are mutually exclusive.
To obtain more detailed information about the movement of the H_{2} molecule inside the cage, the 1D PDF is calculated. Figure 3 shows the 1D PDF, P(R), for (H_{2})_{n}@C_{70} (n = 1, 2) under different temperatures (1 K, 50 K, 100 K, 200 K, and 300 K). Where R is the distance between the center of the hydrogen molecule and the cage center. Figure 3(a) shows that for n = 1, the H_{2} molecule is not fixed in the cage center but has a probabilistic distribution around the cage center. For n = 2, the distance from either of the two H_{2} molecules to the cage center is equal, which may be due to the fact that the two H_{2} molecules are identical and the intermolecular interaction can have an equal effect on them. Figure 3(b) displays the temperaturedependent spatial distributions for two H_{2} molecules in C_{70}, from which we can see when two H_{2} molecules are confined in the cage, their distributions are also probabilistic within a small region. For both n = 1 and n = 2, each P(R) is enlarged and the peak gradually shifts away from the cage center with the increase of temperature. This phenomenon shows that the thermal vibrations of the caged H_{2} molecules are enhanced at a higher temperature. However, in a wide temperature range, R remains much shorter than the average radius of caged C_{70}. Meanwhile, as the temperature increases, the change of P(R) for n = 2 is smaller than that for n = 1, which complementally reflects that (H_{2})_{2}@C_{70} is more stable than H_{2}@C_{70}.
To intuitively compare the spatial distribution differences, the 1D PDF P(R)s for n = 1 and n = 2 under 300 K are plotted together in Fig. 4. For n = 1, H_{2} moves within a region with a radius of 2.5 Bohr and the peak appears near 1.0 Bohr. While for n = 2, the distribution radius of H_{2} molecule varies between 1.0 Bohr and 3.0 Bohr and the peak appears near 2.1 Bohr. Obviously, compared with one H_{2} molecule, two H_{2} have less room to move, so their amplitudes diminish. Furthermore, for n = 2, P(R) = 0 when R < 1.0 Bohr. This is caused by the intermolecular interaction, which stays away from the cage center.
A PIMC method is used to study the thermal properties of H_{2} molecules inside the C_{70} cage. The temperaturedependent interaction energies and spatial distributions of (H_{2})_{n}@C_{70} (n = 1, 2, 3) are calculated. From 1 K to 300 K, the interaction energies are negative for one or two H_{2} molecules caged in C_{70}, while for (H_{2})_{3}@C_{70}, the interaction energy is positive. Therefore, the C_{70} cage can trap two H_{2} molecules at most. Meanwhile the result shows that (H_{2})_{2}@C_{70} is more stable in energy than H_{2}@C_{70}. In both H_{2}@C_{70} and (H_{2})_{2}@C_{70} systems, the spatial distributions show that H_{2} molecule stays further away from the cage center when temperature increases. Meanwhile, this change is slow and the distance from H_{2} molecule to the cage center is much shorter than the average radius of the C_{70} cage. In this paper, for the interaction between H_{2} molecule and the C_{70} cage, we use LJ potential. For further study, a more accurate potential should be used to analyze the rotation of the H_{2} molecules in the cage. Moreover, other properties, such as the spectroscopic properties of (H_{2})_{n}@C_{70} (n = 1, 2, 3) under different ambient conditions, can be studied to find more applications of the systems. Because PIMC is an accurate tool to study quantum manybody properties, it can be used to study other fullerene complexes, such as boron fullerene complexes and boron– nitrogen fullerene complexes.
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 
