Gap opening in the most stable phases of K3 Terphenyl compound

The crystalline structure of potassium intercalated para-terphenyl has been theoretically investigated using a van der Waals Density Functional (vdW–DF) scheme. The incorporation of K atoms into the pristine terphenyl crystal continuously improves the compound stability up to a 3:1 stoichiometry. In addition to confirming previously studied structures, several new crystalline phases showing inedited patterns -small clusters formed by four to six potassium atoms- have been found, in particular, for the more stable K3Terphenyl crystals. Initially, the band structure of these compounds shows metallic features. Nevertheless, a better description of the electronic structure using hybrid functionals that go beyond generalized-gradient-approximation (GGA) reveals that potassium electrons completely filled three sub–bands that develop at the conduction band bottom. Consequently, these stable phases should be discarded as metallic phases able to show superconductivity at low temperatures. Moreover, our study evidences that due to the important electronic correlations shown by these compounds only computational schemes that go beyond GGA are able to provide robust conclusions.


Introduction
The discovery of signs of superconductivity in K 3 Picene (K 3 C 22 H 14 ), with a critical transition temperature (Tc) of up to 18 K in some minor phase [1], has increased the family of high-Tc superconductors to include polycyclic aromatic hydrocarbons (PAH) intercalated with alkali metals (occasionally with other metals), but also has raised some doubts, being received with skepticism from some studies [2]. PAH form molecular crystals with semiconducting or insulating properties which upon thermal intercalation turn into black powders, unstable to air and humidity and so far very difficult to characterize experimentally from a structural point of view. A number of PAH other than picene, including phenanthrene, coronene and dibenzocoronene among others have enlarged the list of materials displaying signs of superconductivity after intercalation in the form of an abrupt diamagnetic change at a given Tc in the magnetic susceptibility vs temperature plot [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The review by Kubozono et al provides an easy starting point to learn about this emerging field [19]. Unfortunately, the alleged superconducting properties rely on a very small fraction of the bulk material, estimated as 1.2% in the case of the K 3 Picene phase above mentioned, and as tiny as 0.04% in the case of K 3 p-Terphenyl [20], object of this work, hindering enormously its proper characterization. Theoretical methods come to help at this point, where experimental data (e.g. X-ray crystallography or Raman spectroscopy among others) give little information about the active phase, largely masked by other overwhelmingly predominant but less interesting phases.
Superconductivity above 120 K was claimed by R-S. Wang and coworkers in March 2017 for a potassium p-terphenyl compound synthesized at a mole ratio of 3:1 [21]. Immediately, experimental evidence for the existence of low energy gaps in a similar material was reported [22]. A few months latter, a possible driving mechanism was published by Mazziotti et al [23]. A detailed study of the magnetization properties of potassium-doped p-terphenyl synthesized at high-pressure was unable to confirm superconductivity despite appearances [13]. Signatures of a possible high-Tc superconducting phase have been also reported in a study of potassium intercalated p-terphenyl flakes obtained by a new procedure [14]. Presently, we are not aware of further support of superconductivity in potassium p-terphenyl compounds.
There are at least two attempts to model the crystal structure of K x Terphenyl using computational methods of Condensed Matter Theory. Both are based on Density Functional Theory (DFT) and they even use the same software package (the Vienna ab initio simulation package, VASP) as we also do. The first one by Zhong et al [24] concludes that formation energies indicate that the system is more stable at x=2.5. On the other hand, the second one by Yan et al [25] founds that the most appropriate structural phase occurs at x=2 and shows a gap of 0.3 eV. Unfortunately, since total energies are not given by them, calculations have to be consistently repeated and formation energies compared. Although this was not the main objective of our work, we have redone geometrical relaxation of both structures with resulting formation energies of −1.75 eV (x=2.5) and −1.67 eV (x=2). Latter, we will show that this is not the end of the story because even better structures can be found for both compositions.
Typically, superconductivity has been claimed for 3:1 compounds. There is a naive reasoning that makes this stoichiometry optimal. It's well known that each potassium atom has one valence electron which can go to the hydrocarbon molecule. Thus, for a 3:1 stoichiometry one has three electrons that will then fill the lowest unoccupied states. The Lowest Unoccupied Molecular Orbital (LUMO) gets completely filled while the LUMO +1 only gets half filled. Therefore, the corresponding band structure would also appear partially occupied pointing to the existence of a Fermi surface and metallic behavior. Many theoretical works based on DFT have shown that this is certainly the case for all studied compounds [26][27][28][29][30][31][32][33]. Nevertheless, the fact that crystalline cells are always occupied by two organic molecules and, consequently, by an even number of potassium atoms for integer stoichiometries suggests that a better electronic arrangement in which electrons are paired would be possible. In this work, we have used hybrid functionals that include a part of the exact Fock exchange to show that this is certainly the case at least for the most stable K 3 Terphenyl structures showing a peculiar kind of potassium clustering. Associated with this surprising result, the charge density characteristics of the occupied sub-bands developing at the conduction band bottom suggests that only two electrons are completely transferred to the organic molecule while the third doping electrons pair two organic molecules thanks to the potassium atoms between them. The impossibility of transferring a third electron to p-Terphenyl is consistent with a deep study recently published [34]. Unfortunately, the use of a non-local potential demands large computational resources and several days or weeks to get converged results. Even so, a wide revision of almost all band structure calculations of metal doped PAHs seems to be necessary.
The rest of the paper is organized as follows. Section 2 is devoted to give some details of the methods and procedures used in this work. Section 3 presents our main computational results together with some discussion of them. The work ends with a few final concluding remarks (section 4).

Computational procedures
All theoretical calculations of crystalline an electronic structures of potassium p-terphenyl compounds been performed using ab initio van der Waals-Density Functional Theory (vdW-DFT) as proposed in 2004 by Dion et al [35]. Taking advantage of an algorithm introduced by Román-Peréz and Soler the whole calculation can be done in reciprocal space [36]. Actually, we have used a later development of the vdW-DFT approximation known as vdW-DFT2 [37]. It was coded by J. Klimeš within the VASP program [38][39][40][41]. Wavefunctions have been expanded in a plane-wave basis set up to a cutoff of 700 eV (that has to be reduced to 400 eV when the exact Fock contribution should be evaluated due to computational limitations) and sampled on a Γ centered Monkhorst-Packgrid automatically generated taking into account the sizes of the reciprocal lattice vectors. One of the files provided within Supplementary Material (KPOINTS) gives the precise recipe consistently used in all computations 3 . The last line of this file is used to obtain the number of points along the reciprocal lattice vectors in such a way that separation between them is about 1/l. We choose l=40 for geometry optimization and l=100 for the density of states calculations resulting in about fifty points in the first case and two to three thousand in the second one. In this way precise densities of states are obtained using the tetrahedra method. Core electrons have been treated within the projector augmented method [42,43]. Semi-core 3s and 3p potassium states have been considered valence states, i.e. they fully enter within the electronic density selfconsistency loop. Consistently, soft potentials have been avoided for hydrogen and carbon, i.e. H, C and K sv have been taken from the POTS file of VASP distribution. Partial occupancies are set by ISMEAR and SIGMA 3 See Supplemental Material atstacks.iop.org/MRX/6/125111/mmedia for a collection of all relevant crystalline structures mentioned in the main text. They are given in Crystallographic Information File (cif) format as defined by International Union of Crystallography. We have systematically used the Structure Data Converter & Editor at Bilbao Crystallographic Server [46] to transform VASP structural data to cif format and vice versa. Additionally, the file KPOINTS generating k-points has been also included. Our inverse cutoff is 40 for geometrical relaxations and 100 for precise DOS calculations. parameters within VASP software. Assuming a metallic behavior at the initial optimization steps, we used ISMEAR=1 and SIGMA=0.2 eV which favors convergence but the scheme was changed to ISMEAR=0 (Gaussian smearing) and very small SIGMA=0.04 eV as soon as a small gap appeared between empty and occupied bands. The relaxation of the electronic degrees of freedom has been stopped when both the total energy and the band structure energy variations between two steps are smaller than 10 −8 eV. Ionic relaxation has been continued as long as any force were larger than 0.001 eV/Å 4 . Unit cells have been always relaxed allowing variations both in shape and volume (ISIF=3 within INCAR file). Volume systematically increases by potassium addition (see table 1) and meanwhile cell shapes completely change as stoichiometry increases. Complete crystallographic information of the studied compounds is provided as Supplementary Material(See footnote 3). The possibility of a breakdown of spin symmetry has been considered in those cases in which very narrow bands are half-filled. This happens seldom since the number of K atoms per cell is even for all integer stoichiometries and also for fractional stoichiometries and doubled or quadrupled unit cells. In some cases that will be discuss later, the possibility of a gap opening appears. In those cases the electronic structure has been calculated beyond GGA using a Hartree-Fock/DFT hybrid functional. Both PBE0 and B3LYP approximations have been employed [44,45]. Notice that in these cases exclusively the electronic states have been recalculated, unit cells and atomic positions were kept invariant. Dispersion terms are not considered within PBE0 or B3LYP functionals since they are based on pure GGA functionals that exclude van der Waals contributions.
We have used free software GRACE maintained by Weizmann to get precise plots representing densities of states [48] and a free version of MERCURY from The Cambridge Crystallographic Data Centre to visualize crystalline structures [49]. Charge isosurfaces have been obtained using XCrySDen software [50]. Structural data from VASP have been translated to Crystallographic Information File (cif ) format using tools provided by Bilbao Crystallographic Server [46].

Results and discussion
Before reporting on details concerning new crystalline phases let us point to some general trends. Usually, large PAH molecules crystallize in phases showing a stacking of layers in which the dominant figure is a herringbone arrangement favoring CH-π binding among organic molecules [51]. Apparently, this kind of molecular arrangement is preserved after the reaction with metals. Experimental structural evidence is scarce [52,53]. Guided by calculated energies of formation of potassium p-terphenyl compounds, we have observed that beyond 2:1 stoichiometry a tendency appears that favors the formation of clusters of potassium atoms formed by four atoms in K 2 Terphenyl, five atoms in K 2.5 Terphenyl, and six atoms in K 3 Terphenyl. These structures are based on a doubled cell in the layer planes formed by organic molecules. This increased unit cell contains four p-terphenyl molecules. Figures 1, 2 and 3 show these kind of layers for increasing potassium content. A similar structure for Cs 2 Phenanthrene is shown and discussed in figure 1 of the work by Takabayashi et al [52]. Table 1. Main properties of the some of the possible compounds formed by the intercalation of potassium into the structure of p-terphenyl polycyclic aromatic hydrocarbon. Only a small fraction of the studied structures has been included in this table (See footnote 6). Two structures previously reported as the most stable ones have been re-optimized using our own recipe to allow a detailed comparison. Some structural details, cell volume and total energy per cell are given for each compound ending with the normalized energy of formation corresponding to two p-terphenyl molecules and the appropriate number of potassium atoms.  4 This set of parameters together with the tag accurate within the VASP input seems to work correctly even for compounds formed by a mixture of organic molecules and alkaline elements. Actually, the selected energy cutoff of 700 eV is well above the recommended values for C, H, Cs and K PAW potentials which are 400, 250, 220 and 259 eV, respectively. Two realistic tests discussed in [47], polymorphism in pristine pentacene and the structure of the KC 8 intercalation compound, were shown to give satisfactory result using the same numerical procedure.
The starting seeds for a compound of a given composition are mainly based in already optimized structures both containing less or more potassium. Once, the experimental structure of K 2 Picene [53] was adapted to K 2 Terphenyl simply removing some carbon atoms and adjusting hydrogen saturation. Afterward it was ionically  relaxed. In other cases, some structures based on chemical reasons were exceptionally used as starting points for new structures. In these way, about one hundred possibilities have been explored. It goes without saying that this procedure does not exhaust in any way all possible structures. Beside foreseeable or already studied patterns, several new phases have been discovered. In particular, all the crystals showing the most stable structures have not been previously reported, except for the one corresponding to KTerphenyl alloy [24] 5 . Formation energies of the studied compounds are obtained from total energy values by a simple difference: where always the original pristine cell of terphenyl formed by two organic molecules is the reference. Therefore, in those cases in which the unit cell of the new stable structure contains four or even eight organic molecules, total energies are normalized to values corresponding to two terphenyl molecules, that is, they are divided by two or four, respectively. Notice also that total energies as given by VASP code are internal energies in the zero temperature limit, volume changes do not contribute to them as they would to enthalpies 6 . Table 1 summarizes our results for the structures showing better stability for each alloy composition 7 . Although formation energies are the relevant physical magnitudes characterizing compound stabilities, total energies obtained with a reproducible suite of calculations (VASP in our case) should help researchers to make precise comparisons between structures obtained by different groups. Actually, these data are usually given in Quantum Chemistry papers.
Our computations have proceeded in two stages. In the first stage, simulations are based on a conventional GGA approximation supplemented by dispersion corrections as given by the van der Waals Density Functional of Langreth and Lundquist groups [35,37]. Actually, we use vdW-DF2 version that was implemented in VASP by Klimeš et al [38] based on the algorithm by Román-Pérez and Soler [36]. This is a well-known path that allows relatively quick simulations that can search for all imaginable structures. In this way, good structures have been obtained for all stoichiometries. The second stage is reserved for those cases in which relevant changes in  6 Internal energies can be safely used in solid reactions instead of the more rigorous enthalpy of formation because volume changes are minimal. For example, the enthalpy variation associated to the volume change occurring when K 3 Terphenyl is formed (penultimate row of table 1) amounts only about 10 −4 eV. 7 We leave for a future publication a comprehensive report on the more relevant structures of K x Terphenyl that we have found. Some of them could be decisive for the explanation of superconductivity based on pseudo-stable structures or for the explanation of x-ray experiments based on potassium quasi-segregation. the calculated density of states can be anticipated. This happens when potassium clustering leads to clusters formed by six atoms and the DOS outlines the formation of sub-bands at the bottom of the conduction band. In such conditions, the effect of the inclusion of a part of the exact exchange energy has been considered. To this end, hybrid Hartree-Fock/DFT functionals have been used. Precise calculations of this type need between 10 and 100 times more computational resources than conventional DFT calculations. This is the reason for its low implantation in Materials Research. Actually, searching for new structures using this formalism is prohibitive these days. First stage results for compounds of increasing stoichiometry follow.

KTerphenyl
Two potassium atoms enter in the terphenyl unit cell without changing the herringbone framework. Nevertheless, space symmetry is reduced because the potassium atoms that enter in nearest-neighbors channels occupy non-equivalent positions. This structure was obtained in 2018 by Zhong et al [24]. The complete symmetry of the pristine terphenyl crystal can be recovered at an energetic cost of 0.37 eV per cell. Let us comment that smaller contents of potassium can also be shown to form stable structures but their interest is small. For example, K 0.5 Terphenyl has potassium occupying alternating herringbone channels, slightly distorted unit cells and a formation energy of −0.3 eV per cell.

K 1.5 Terphenyl
It could be expected that the third potassium atom entering the unit cell would add to one of the two herringbone channels occupied by one alkaline atom. This is certainly possible but not the best option from an energetic point of view. Just for this composition, the compound shows the tendency that is crucial in our work which is none other than the tendency of potassium to form atomic aggregates. Three K atoms fully occupy one of the channels leaving free the second one. The k1.5terphenyl.cif file supplied within Supplementary Material contains complete crystallographic information(See footnote 3). Certainly, a more regular distribution of potassium in groups of one and two atoms is also possible but has an energetic cost of 0.1 eV. By the way, this 2-1 pattern is the one given in [24].

K 2 Terphenyl
The aggregation tendency of potassium atoms shows once again at this 2:1 stoichiometry. Instead of entering two atoms in each herringbone channel as it has been found for other alloys of the same kind 8 , the crystalline structure evolves towards a doubled cell herringbone structure with four terphenyl molecules per cell and potassium clusters formed by four atoms. As it happened for K 1.5 Terphenyl, only half of the new herringbone channels appear occupied as figure 1 shows. A complete knowledge of the whole crystalline structure can be obtained from k2terphenyl.cif file given in Supplemental material(See footnote 3). Terphenyl molecules show large distortions when compared with free copies. Let us notice that the formation energy of the more conventional herringbone structure is −1.565 eV, that is, only 0.138 eV less stable per unit cell. Important distortion of terphenyl molecules are absent in this well studied structure. Let us also comment that many more investigated structures give competing energies of formation for this particular composition including not so regular herringbone structures with three and one K atoms in the channels and also slipped parallel patterns [47,54]. In any case, corresponding DOS show completely occupied bands, that is, a gap separating the occupied part of the conduction band from the empty part. Compounds are predicted to be insulators.

K 2.5 Terphenyl
The most stable structure for a 5:2 stoichiometry shows an even larger unit cell that contains now eight terphenyl molecules. Although the overall arrangement resembles the one shown in figure 1, the number of K atoms populating the channels varies here from four to six (see figure 2). It is actually the only way to reach this Figure 5. Density of states (DOS) of pristine terphenyl (a) and 3:1 potassium terphenyl compound in its most stable structure as obtained for four different theoretical models: (b) usual vdW-DF2 results, (c) pure GGA results using PBE recipe, (d) PBE0 results obtained for an hybrid functional that includes a 25% exact Fock exchange, and, (e) B3LYP results obtained using one of the most popular hybrid functional of Quantum Chemistry (it includes a 20% exact exchange). Fermi levels are indicated by black dashed lines. A gap separating unoccupied levels from occupied ones appears as soon as exact exchange is considered. Each of the sub-bands at the beginning of the conduction band holds four electrons, twelve in total as the number of K atoms that are present in the double herringbone cell. particular mixing ratio. In this case, the best formation energy is 0.22 eV per cell more bonding that the one obtained in [24] using a conventional unit cell containing two terphenyl molecules.

K 3 Terphenyl
Mixing of potassium and p-terphenyl reaches its maximum energy of formation precisely at this 3:1 stoichiometry with a value of 2.2-2.3 eV per unit cell depending on details (see table 1). This is about 0.3 eV better than K 2.5 Terphenyl and more than 1 eV when compared with K 4 Terphenyl (not included in the mentioned table) 9 . Figures 3 and 4 show the two found arrangements of planes in this case. Here p-terphenyl molecules show less distortion relative to the free case than the one obtained for lower stoichiometries. The main feature present in both crystals is the formation of well-ordered potassium clusters containing six atoms. The number six is very relevant since previous experimental results obtained by Carrera et al [34] suggest that no more than two electrons can be transferred from K atoms to p-terphenyl molecules. If these is the case, a cluster formed by six K atoms would transfer four electrons to the aromatic molecules (two electrons per molecule) remaining two electrons in the neighborhood of the cluster in a deep state. The important point is that two electrons allow pairing which is impossible in alternative lattices like pristine-based crystals as shown in figure 4(e) of [24]. This novelty opens the possibility of fully occupied bands even for this 3:1 stoichiometry. Therefore, we have calculated the Density of States (DOS) in both cases. Results are shown by figures 5(b) and 6(a). Both DOS show the absence of gaps within the conduction band. A metallic behavior would be predicted in line with already published work. We leave to the next subsection the improvement of this calculation to better include correlation effects.

K 6 Terphenyl
We have obtained a well converged crystalline structure for this compound based on the x=3 results. Here large clusters containing twelve potassium atoms are formed. They occupy the same channels shown in figure 3 Here each band holds two electrons. 9 We have work hard to find the best structure for K 4 Terphenyl compound without fully achieving the proposed goal. In any case, our still unfinished results clearly show that x=3 provides the largest stabilization energy. but lattice constants are considerably larger. Detailed information is given by k6terphenyl.cif file in Supplementary Material(See footnote 3). The formation energy of this compound confirms both the saturation of intercalation beyond x=3 and the tendency of potassium to form clusters.

3.2.
Beyond GGA: Hybrid Hartree-Fock/DFT functionals Second stage calculations have been restricted to K 3 Terphenyl structures because numerical simulation shows that this is the most favorable composition and DOS predicts metallic behavior. In principle, they are candidates for superconducting phases. On the other hand, optimized K 2 Terphenyl structures show a gap [25] which would probably increase with a better description of electronic correlation. K 2.5 structures are difficult challenges. Both are metallic at the GGA level and, therefore, candidates for high-Tc superconductivity [24]. Our cell contains eight terphenyl molecules which make the use of hybrid functionals impossible (see figure 2). On the other hand, the cell of the not so stable structure shows a five atoms cluster where the atomic positions of K atoms are not completely fixed. We guess that disorder could play an important role in this case. If this were the case, results obtained for a static configuration would probably be not relevant. For all these reasons, only K 3 Terphenyl structures are further analyzed.
Figures 5(b) and 6(a) show DOS predicting a metallic behavior of the corresponding 3:1 compounds. Actually, this would be a good starting point for the appearance of superconductivity at this stoichiometry. Moreover, the situation coincides with a naive filling of the molecular orbitals of p-terphenyl molecules: LUMO is completely filled with two electrons but LUMO+1 is half occupied. Nevertheless, the existence of clusters formed by six K atoms could suggest that two electrons would remain trapped in the cluster at a new mixed potassium-terphenyl deep orbital. If these were the case, an additional band could detach from the lower part of the conduction band. In order to explore this possibility a more accurate description of correlation among electrons is necessary. In consequence, DOS have been recalculated including part of the exact exchange. One should abandon van der Waals contributions included in vdW-DF2 and start from a pure GGA functional. Inside VASP, PBE is the natural choice. Figures 5(c) and 6(b) show our GGA results. Interestingly, three incipient sub-bands appear at the lower part of the conduction band although no gaps are still formed. The next step is the use of some hybrid Hartree-Fock/DFT functional. Again, within VASP the natural choice is PBE0 that includes a 25% exact Fock exchange. Computing becomes very demanding for these non-local calculations. Runs are typically between 10 an 100 times slower even if good parallel resources can be accessed. In this situation it is not possible the refinement of the atomic structure because just the computation of the electronic structure for a static atomic configuration takes a long time. Results for PBE0 model are shown in figures 5(d) and 6(c). Both show three well separated bands at the bottom of the conduction band. More importantly, a gap divides the occupied part of the DOS from the empty one pointing to an insulating behavior at low temperatures. Gaps extend over 0.27 and 0.14 eV for the structures shown in figures 3 and 4, respectively. In the first case, we have been able to check that the particular model used to include exact exchange does not change the new paradigm. Using b3lyp, a very popular hybrid functional in Quantum Chemistry studies, we get a completely consistent result as figure 5(e) shows.
The analysis of the wavefunctions involved in the formation of the conduction band bottom confirms our initial guess. Charge isosurfaces have been separately obtained for the three sub-bands at the lower part of the Figure 9. Charge isosurfaces obtained after subtracting the results obtained for artificial crystals formed exclusively by either potassium atoms or terphenyl molecules in the same arrangement as the compound. Comparison with a typical K-K isosurface. conduction band (see figure 7). With some caution the first and second bands can be assigned to LUMO orbitals (left panel of figure 8) of p-terphenyl molecules while the third one shows a more complex aspect that we ascribe to some resonance between LUMO+1 organic orbitals (right panel of figure 8) and 4s K orbitals. When the charge distribution of virtual crystals containing exclusively either terphenyl molecules or potassium atoms is subtracted from the charge distribution of the real system, the result shown in figure 9 is obtained. One can see that charge flows from K neighborhood to the terphenyl region. The inset shows a typical K-K charge distribution at the equilibrium distance, that is, for much closer atoms. Based on this analysis, we can conclude that although potassium only transfers two electrons per organic molecule, the existence of six-atoms clusters provides a good pairing possibility for the two remaining electrons.

Concluding Remarks
Lets recapitulate the main findings of our work: (i) Potassium terphenyl compounds showing 3:1 stoichiometry are more stable than 5:2 compounds. Its crystalline structures are characterized by the formation of clusters containing six potassium atoms. (ii) Once electronic correlation is described beyond GGA, a gap opens within the lower part of the conduction band that separates occupied states from empty ones. (iii) Consequently, an insulating low temperature behavior is predicted for the more stable 3:1 compounds. (iv) A prediction of metallic behavior is restricted either to the not so stable K 2.5 Terphenyl compound showing clusters of five potassium atoms or to small deviations from exact 3:1 stoichiometry. In both cases, we believe that the intrinsic disorder inherent to this cases would make predicted metallicity very doubtful. (v) We think that the study of fluctuations of potassium 3s and 3p deep levels using XPS could help to the knowledge of the actual existence or not of potassium clustering in these materials.