First-principles study of the charge ordered phase in $\kappa$-D$_3$(Cat-EDT-TTF/ST)$_2$: Stability of $\pi$-electron deuterium coupled ordering in hydrogen-bonded molecular conductors

TL;DR

This study uses first-principles DFT calculations to investigate the charge-ordered phase in hydrogen-bonded molecular conductors, revealing the importance of hybrid functionals for accurate structural predictions and identifying a stable noncentrosymmetric charge-ordered phase.

Contribution

It demonstrates that hybrid functionals are essential for accurately modeling charge-ordered phases in hydrogen-bonded molecular conductors, and uncovers a new stable noncentrosymmetric charge-ordered structure.

Findings

Hybrid functionals reproduce charge-ordered phase stability.

Standard DFT fails to capture structural stability.

A new noncentrosymmetric charge-ordered phase is identified.

Abstract

We study the electronic and structural properties of the low-temperature ordered phase of hydrogen-bonded molecular conductors, $\kappa$-D$_3$(Cat-EDT-TTF)$_2$ and its selenium-substituted analog $\kappa$-D$_3$(Cat-EDT-ST)$_2$, by means of first-principles density functional theory~(DFT) calculations. In these compounds, the charge ordering in the $\pi$-electron system is coupled with the ordering of the displacements in the deuteriums forming the hydrogen-bond, equally shared by two oxygens in the high-temperature phase. While the structural optimization within the standard DFT method based on the generalized gradient approximation fails to reproduce the structural stability of the charge-ordered (CO) phase, we show that a hybrid functional of Heyd, Scuseria, and Ernzerhof can reproduce structural characters of the CO phase, owing to the more localized nature of the wave functions. Furthermore, using the ability of the hybrid functional to predict the electronic and structural properties, we find a stable noncentrosymmetric CO phase with another pattern of deuterium ordering.

Tables

Table 1. Table 1: Experimental and theoretically optimized structural parameters of central C=C bonds in the TTF part and hydrogen bonds in κ 𝜅 \kappa -D 3 (Cat-EDT-TTF) 2 (D-S) and κ 𝜅 \kappa -D 3 (Cat-EDT-ST) 2 (D-Se). The P 𝑃 P 1 structure for D-S is the noncentrosymmetric structure discussed in the text.

	Method	C=C [Å]	C=C [Å]	O$\cdots$O [Å]	O–D [Å]	O–D [Å]	$\angle$O–D–O [deg]
		w-D	w/o-D		w-D	w/o-D
D-S
Centrosymmetric ($P\overline{1}$)	Exp.	1.34	1.40	2.50	1.02	1.51	164.7
	GGA-PBE	1.37	1.37	2.43	1.14	1.29	173.9
	HSE06	1.35	1.38	2.48	1.06	1.41	173.9
Noncentrosymmetric ($P1$)	HSE06	1.37	1.38	2.43	1.10	1.33	174.7
		1.35	1.36	2.49	1.05	1.43	173.3
D-Se
Centrosymmetric ($P\overline{1}$)	Exp.Ueda et al. (2015)	1.35	1.38	2.52	0.80	1.73	167.9
	GGA-PBE	1.37	1.37	2.43	1.09	1.40	172.9
	HSE06	1.35	1.38	2.48	1.06	1.43	173.8

Full text

††preprint: APS/123-QED

I Introduction

Accurate calculation of structural and electronic properties for crystalline solids hosting strongly correlated electron systems is a long-standing problem in condensed matter physics. The standard first-principles calculations based on the density functional theory (DFT) which have shown great success in many materials often fail to reproduce the insulating states such as Mott insulatorsTerakura et al. (1984); Miyazaki and Ohno (1999); Nakamura et al. (2009) and charge ordered (CO) insulators. Yamauchi et al. (2009); Yoshimi et al. (2012); Jacko et al. (2013); Yamakawa et al. (2016); Katoh et al. (1999); Anisimov et al. (1996); Jeng et al. (2004); Botana et al. (2016) One of the prototypical systems is molecular conductors in which correlated electrons are formed owing to comparable energy scales for the kinetic energy and the Coulomb interactions among electrons.Seo et al. (2004) A fundamental problem is that the CO phases in these molecular systems, often seen experimentally and extensively studied theoretically, are not structurally stable within the standard DFT approach. Namely, even when we perform structural optimization starting from the experimentally observed low-temperature structure with charge ordering, the optimized structure becomes that of the high-temperature phase where the charge disproportionation is absent.

Here we study a class of molecular crystals based on catechol with ethylenedithiote-tetrathiafulvalene, Cat-EDT-TTF, and its selenium-substituted analog, Cat-EDT-ST. Kamo et al. (2012) These systems are characteristic in the sense that their electronic and structural properties are strongly linked, through the hydrogen bonding in their constituent molecular units. The units form the so-called $\kappa$-type two-dimensional arrangement; the solids are called $\kappa$-H3(Cat-EDT-TTF)2 (hereafter abbreviated as H-S) and $\kappa$-H3(Cat-EDT-ST)2 (H-Se). Their deuterated samples, $\kappa$-D3(Cat-EDT-TTF)2 (D-S) and $\kappa$-D3(Cat-EDT-ST)2 (D-Se), are also prepared and experimentally investigated. Interestingly enough, large differences between the hydrogen and deuterium samples are seen. At high temperatures, they all show Mott insulating behavior with localized $S$ = 1/2 spins. Isono et al. (2013, 2014) Structurally, every two H/D(Cat-EDT-TTF/ST) units share a H/D atom forming the hydrogen-bonding between two nearest oxygens, bridging the molecules with a relatively short O$\cdots$O distance [Fig. 1(a)]. These are consistent with the electronic structure calculations Isono et al. (2013); Tsumuraya et al. (2015) where the $\pi$-electrons form a half-filled band owing to the dimerization in the $\kappa$-type arrangement of H/D(Cat-EDT-TTF/ST) units which are all equivalent. The strong electron-electron repulsion is the probable reason the system is Mott insulating (the so-called dimer-Mott insulatorSeo et al. (2004); Kino and Fukuyama (1996); Miyazaki and Ohno (1999); Nakamura et al. (2009)). Then the $S$ = 1/2 spins are localized on every dimer, namely the antibonding pair of the highest occupied molecular orbitals (HOMO) of H/D(Cat-EDT-TTF/ST) molecules.

The difference between the H and D samples appears at low temperatures. The H samples stay paramagnetic down to the lowest temperatures, showing the possibility of a quantum spin liquid state. Isono et al. (2014) The coupling between the quantum fluctuation of protons and $\pi$-electrons is discussed based on a dielectric constant measurement showing a quantum paraelectric behavior. Shimozawa et al. (2017) In fact, our previous first-principles DFT calculations show that the optimized distance between the two O atoms is relatively short, and the calculated potential energy surface for the shared H atom is very shallow near the minimum points. Tsumuraya et al. (2015)

In stark contrast, in the D samples, a first-order phase transition occurs that is associated with a large structural change at a transition temperature ($T_{c}$) of 182 and 185 K in D-S and D-Se, respectively. Ueda et al. (2014, 2015, 2019) The magnetic susceptibility decreases sharply below $T_{c}$ resulting in a nonmagnetic ground state. Shimozawa et al. (2017); Yamashita et al. (2017) The heat capacity of D-S shows smaller values than in H-S at low temperatures, consistent with the absence of spin contribution.Yamashita et al. (2017) Below $T_{c}$, the D atom forming the hydrogen bond localizes near one of the two O atoms. The low-temperature structure (space group: $P\bar{1}$) suggests the existence of charge disproportionation between two types of dimers that consist of two distinct molecules, D2(Cat-EDT-TTF/ST) and D(Cat-EDT-TTF/ST), abbreviated here as w-D and w/o-D units, respectively [Figs. 1(b) and (c)]. Ueda et al. (2014, 2015) It is considered that the charge ordering in the $\pi$-electron system is coupled with the D ordering. Naka and Ishihara (2018)

In our previous study, Tsumuraya et al. (2015) we performed structure optimization for H-S and found a H-coupled CO phase, which has the same symmetry as the low-temperature structure of D-S and D-Se. However, the results show a metallic state with a large Fermi surface in contrast to the insulating behavior in experiments, and the magnitude of charge imbalance, and the structural distortions from the high-temperature phase are much smaller than those observed in the D samples. These calculations were done with a conventional exchange-correlation functional of the generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof (PBE), Perdew et al. (1996) which often underestimates the localization of electrons and then the structural stability of CO phases. Since the band structure, distortion of molecules, and their packing are all sensitive to the localization of electrons and the magnitude of the charge imbalance, GGA may not have enough accuracy to calculate the electronic structures of such molecular CO systems. Due to this problem of GGA-PBE, it is very difficult to construct a reliable effective model, which is essential for many theoretical studies related to the CO state. Some of the parameters in the effective model may be determined from the experimental results, but information from first-principles calculations is highly desirable.

In this work, we compare the GGA-PBE functional with another exchange-correlation functional, i.e., a hybrid functional by Heyd, Scuseria, and Ernzerhof (HSE06), Heyd et al. (2003); Heyd and Scuseria (2004); Heyd et al. (2005, 2006) and investigate the stability of the CO states in D-S and D-Se. We expect that the hybrid functional method will provide more reliable results than GGA to calculate the electronic and structural properties of this class of materials. Giovannetti et al. (2009); Fonari et al. (2014) However, since the computational cost is much more expensive than that of GGA, they have not been applied to complex molecular CO systems and their quantitative accuracy is yet to be evaluated.

The rest of the paper is organized as follows. We introduce the crystal structure and the calculation method in Secs. II and III, respectively. In Sec. IV, we show the results: The electronic structure is discussed based on the experimental structures of D-S and D-Se and the use of GGA-PBE and HSE06 functionals are compared (Sec. IV A). Structural optimization for D-S is performed, and the structural stability of the CO phase is discussed, including the possibility of another CO pattern (Sec. IV B). Sections. V and VI are devoted to discussion and a summary, respectively.

II Crystal structure

The low-temperature phase of D-S and D-Se has a triclinic structure (centrosymmetric space group: $P\bar{1}$) where one of the D atoms is located at an off-center position closer to one of the two O atoms forming hydrogen bonds, as shown in Fig. 1(a). This is in contrast to the high-temperature phase (centrosymmetric space group: $C2/c$), isostructural to the H samples, where all the H/D(Cat-EDT-TTF/ST) units are equivalent and the hydrogen bond-forming H/D is equally shared by two oxygens in different units. Figures 1(b) and 1(c) show the crystal structure, which is characterized by a charge disproportionation between monomers w-D and w/o-D. Each monomer forms a dimer with the same type of monomer; namely, there are w-D dimers and w/o-D dimers with a 1:1 ratio. These two kinds of dimers, which are in the high-temperature $C$2/$c$ structure equivalent and connected by glide operations, are alternately stacked with a periodic arrangement.

III Calculation method

In order to study the electronic and structural properties of the CO state coupled with D ordering, we compare the results with two first-principles approaches mentioned above, using GGA-PBE and HSE06 as exchange-correlation functionals. Note that we need to treat strong intramolecular and weak intermolecular interactions simultaneously. Furthermore, the energy difference between competing phases is often very small, and can be as small as of the order of 10 meV/f.u. in our case, as discussed later.

For the GGA calculations, Kohn-Sham equations are self-consistently solved in a scalar-relativistic fashion using the all-electron full-potential linearized augmented plane wave (FLAPW) method implemented in the QMD-FLAPW12 code Wimmer et al. (1981); Koellng and Arbman (1975); Weinert (1981); Tsumuraya et al. (2012). The band structure calculations with GGA are also performed using the pseudopotential method based on the projector augmented wave (PAW) method with plane wave basis sets implemented in the Vienna $Ab$ $initio$ Simulation Package. Vanderbilt (1990); Kresse and Furthmüller (1996); Kresse and Joubert (1999) The results are fairly in agreement with each other. The LAPW basis functions in the interstitial have a cutoff energy of 30.3 Ry. The angular momentum expansion inside the sphere is truncated at $l$ = 8 for all the atoms. The cut-off energy for the potential and density is 276 Ry. The muffin-tin sphere radii are set as 0.66, 0.36, 0.62, 1.03 and 1.16 Å for C, H, O, S, and Se atoms, respectively. The $\bf{k}$-point meshes used are 4 $\times$ 4 $\times$ 2.

As for the HSE06 hybrid functional, the Kohn-Sham equations are solved using the pseudopotential method based on PAW method with plane wave basis sets implemented in VASP. Paier et al. (2006) Previously, the role of the exact exchange on the electronic structure for the charge modulated state was discussed for TiSe2 Chen et al. (2015); Hellgren et al. (2017) and BaBiO3, Foyevtsova et al. (2015); Smolyanyuk et al. (2017) and it was shown that the use of the hybrid functional was essential for obtaining a proper description of the electronic and structural properties. In the present HSE06 calculations, we first obtain a converged charge density from the self-consistent calculation within GGA, and then the self-consistent hybrid functional calculations are performed using the GGA charge density as the initial state. A common $\bf{k}$-point sampling, including for the structural optimization, is set as 3 $\times$ 4 $\times$ 2. The cut-off energy for plane waves is 29.4 Ry for the HSE06 calculations. The range-separation parameter in HSE06 calculations is 0.2 Å*-1*, and 25$\%$ of the exact exchange is mixed to the GGA exchange for the short-range interactions.

We note that the present study does not consider the quantum effects of H or D atoms. Ishii et al. (2006); Horiuchi et al. (2010); Ishibashi et al. (2018) We do not distinguish H and D atoms in our DFT calculations, and the H/D isotope effect is not discussed.

For the lattice parameters, we use the experimental ones throughout this paper, since we have difficulty with their theoretical determination as explained in the following. First, the importance of dispersion interactions or van dar Waals interactions, which are not correctly treated in the GGA-PBE or HSE06 functional, in the CO phase is unclear. D3_ We expect the electrostatic interactions in the CO phase are much stronger. Seo et al. (2004); Naka and Ishihara (2018)

Second, if the effect of dispersion interactions is important, the choice of the method treating the dispersion corrections would significantly affect the optimized lattice parameters. To evaluate the lattice parameters accurately, we need to use a method which can treat the dispersion interactions as well as the CO state, simultaneously. We may also have a problem in the numerical accuracy to calculate the stress tensors since they are subtle quantities related to the calculation conditions. Considering these ambiguities in the optimization of lattice parameters and, in contrast, the high accuracy of the experimental ones especially at ambient pressure, we use the experimental lattice parameters and concentrate on clarifying the detailed accuracy of the first-principles evaluation of various quantities related to the CO state of the title compounds. Since the reliability and accuracy of first-principles methods within GGA-PBE or HSE06 for the CO state in molecular conductors are far from clear, we believe that our present study is also valuable for constructing a reliable effective model.

IV Results

As explained above, it is important to investigate the accuracy of the first-principles methods for the title compounds. In this section, we first investigate the electronic properties of the CO state using the experimental structure, including the internal coordinates, in Sec. IV A. Then, in Sec. IV B, we work on its structural stability by optimizing the internal coordinates using GGA-PBE and HSE06 functionals.

IV.1 Electronic structure

In this section, the electronic structure based on the experimental crystal structure in the CO state is discussed. Figure 2(a) shows the calculated band structure of D-S within GGA. The top of the valence band is located at the $Y$-point and the bottom of the conduction band is located at the $U$-point, which both slightly cross the Fermi level; the system is (semi-)metallic.

As plotted in Fig. 2(b), the calculated local density of states (LDOS) clearly shows charge disproportionation between the two types of monomers in the unit cell. The solid and dashed curves are the LDOS of the w-D and w/o-D units, respectively. The LDOS are obtained as a summation of projected DOS on C $p$, O $p$, S $p$ and S $d$ states in the respective monomer units. Compared with the CO state discussed in our previous study for H-S, Tsumuraya et al. (2015) the magnitude of the intermolecular charge disproportionation is much larger, which is consistent with the experimental observations. At around -0.1 to 0 eV and –0.5 to –0.3 eV, where the second and third bands locate among the four bands, the LDOS of w-D is much larger than that of w/o-D. Conversely, the top and fourth bands at around 0 to +0.2 eV and –0.7 to –0.5 eV, respectively, mostly originate from the LDOS of w/o-D.

These features can be understood by a schematic energy diagram of molecular orbitals of the two kinds of dimers, as shown in Fig. 3. Here, the energy levels of monomers are evaluated by the calculations of the isolated monomers in a supercell. We first note that the HOMO level of the w/D unit is lower than that of the w/o-D unit. This is reasonable since we expect that it is stabilized by the hydroxyl (-OD) group at the end of the molecule.

When they form dimers, the energy levels are split into antibonding and bonding states, and since the distances between the two monomers are very different [see Figs. 1(d) and 1(e)], the splitting is much larger in the w/o-D units with shorter distance (3.42 Å) than in the w-D units with larger distance (3.71 Å). These MO levels of the dimers are obtained from a single-$\bf{k}$-point calculation for the crystalline solid. As a result, the antibonding level of the w/o-D dimer is higher than that of the w-D dimer, and the occupancies of the two antibonding levels become different, 2 for w-D and 0 for w/o-D dimers, consistent with the realization of the CO state in the solid. From these analyses, we conclude that there are two main contributions to stabilize the CO state: i) the energy difference between w/o-D and w-D monomers, and ii) the stronger dimer formation of w/o-D than w-D. Since the former is closely related to the localization of H or D atoms, their positions and dynamics are coupled to the stability of the CO state.

Next, the GGA band structure of D-Se and the LDOS of the w-D and w/o-D units are shown in Figs. 2(c) and 2(d), respectively. Similarly to D-S, the Fermi level crosses the upper two bands resulting in a metallic state, and the LDOS shows charge disproportionation between the w-D and w/o-D units. The difference from D-S is that the dispersion of each band is about 15 $\%$ larger than that of D-S while, from the structural point of view, the interplanar distances between monomers are slightly longer [Figs. 1(d) and  1(e)]. These originate from the fact that the Se-$p$ state is more delocalized than S-$p$. Again, good correspondence between the four bands and the LDOS of the four monomer units is seen. Therefore, we expect that the mechanism of CO formation coupled to the D ordering is the same as in D-S.

Now, the band structures and the LDOS of D-S and D-Se using the HSE06 hybrid functional are plotted in Fig. 4. Compared to the GGA results, the four bands are farther apart from each other, suggesting the more localized behavior of the wave functions, as expected. Since the same correspondence between the bands and the molecular orbital levels discussed above holds, the enhanced separation of the bands leads to more pronounced charge disproportionation. This is seen in the LDOS in Figs. 4(b) and 4(d), especially in the upper two bands where the LDOS of the antibonding state of w-D (w/o-D) units becomes more occupied (unoccupied) compared to the GGA case in Fig. 2.

In D-S, more importantly, a finite band gap is opened. An indirect band gap of 0.04 eV is obtained between the top of the valence bands at the $Y$-point and the bottom of conduction bands at the $U$-point. The same tendency is seen in the band structure of D-Se, but it is not enough to make the system an insulator. Although there are some quantitative differences, GGA and HSE06 provide similar band structures for both compounds. However, we will show that the structural stability of the CO state is completely different between these two functionals, in the next section.

IV.2 Structural optimization and possible noncentrosymmetric phase

In this section, we investigate the accuracy of structure determination of the internal coordinates using the experimental lattice parameters, as in our previous work. Tsumuraya et al. (2015) First, we discuss the results obtained by setting the initial state as the experimental structures of the D localized phase in D-S and D-Se. We have also found another stable structure with noncentrosymmetric space group $P$1, which will be discussed later.

The experimental and theoretically optimized bond parameters are summarized in Table 1, and the O$\cdots$O and C=C distances are shown in Fig. 5. The first thing we notice is that, after the geometrical relaxations within GGA-PBE, the difference in the C=C bond length at the center of the TTF part of w-D and w/o-D is less than 10*-3*Å, in both D-S and D-Se. This is also the case for the rest of the bond lengths in the two monomer units, suggesting that the structure is basically relaxed to that of higher symmetry $C$2/$c$ of the high-temperature phase, within the numerical accuracy. Since this C=C bond length is known to be sensitive to the charge disproportionation, the results indicate that CO state is unstable after the structural optimization.

We can overcome this problem using the HSE06 hybrid functional. Now, as seen in Table 1 and Fig. 5, the optimized structures are fairly in agreement with the experiments. Importantly, the central C=C bond lengths in the TTF part are well reproduced, 1.35 and 1.38 Å for the w-D and w/o-D units (for both D-S and D-Se) respectively, within $\pm$ 0.02 Å compared with the experimental values. Therefore, we conclude that the HSE06 functional is highly accurate for reproducing the structures of the molecular CO phase quantitatively.

Encouraged by the reliability of the results above, we now investigate another CO state with D ordering, but with a different pattern from the experimental observations. We consider the CO pattern known as the “ferroelectric” phase discussed in $\kappa$-(BEDT-TTF)2$X$ BED and related systems. Naka and Ishihara (2010); Abdel-Jawad et al. (2010); Gomi et al. (2010); Naka and Ishihara (2013); Gomi et al. (2013); Gati et al. (2018); Hassan et al. (2018) It gives rise to a noncentrosymmetric $P$1 structure, and its possible realization coupled to the displacement of H/D in H/D-S has been discussed in Ref. Naka and Ishihara (2018) based on an effective model.

Here, to model the noncentrosymmetric phase, in the initial state we change the position of D atoms, from the experimental structure of the low-temperature phase of D samples, to that where the direction of the displacement of D from the middle point between the shared oxygens is opposite. We indeed find a stable solution by structural optimization, as listed in Table 1. One can see that there are four kinds of C=C bond lengths as well as O-D distances, which we can classify into two each: the w-D (w/o-D) units with shorter (longer) C=C bonds, indicating larger (smaller) electron occupation. The pattern is shown in Fig. 6, together with the experimentally observed pattern [the same as in Figs. 1(b) and (c) but with expanded views]. In the latter [Figs. 6(a) and 6(b)], the two kinds of dimers, i.e., the w-D dimers and w/o-D dimers are arranged in a checkerboard pattern in the two-dimensional planes. On the other hand, as shown in Figs. 6(c) and (d), in the noncentrosymmetric phase, each dimer has one w-D unit and one w/o-D unit, leading to charge disproportionation and electric polarization within the dimers. The pattern globally breaks the inversion symmetry.

Energetically,  the centrosymmetric and noncentrosymmetric structures are very close, while the centrosymmetric phase (-406.190 eV/f.u.) is slightly more favorable than that of the noncentrosymmetric structure (-406.182 eV/f.u.) within an energy difference of 8 meV per formula unit.  It is an interesting issue to evaluate the value of electrical polarization and identify its nature in the noncentrosymmetric phase Horiuchi et al. (2010); King-Smith and Vanderbilt (1993); Resta (1994); Ishibashi et al. (2010), under the variation from a “paraelectric” phase involving both D displacements and molecular deformation, which we leave as a future problem.

V Discussion

Using the HSE06 functional, our theoretical results by structural optimization show overall good agreement with the experimental structural parameters of the low-temperature phase in D-S and D-Se. This indicates that the more localized nature of the resulting wave functions is reflected in the electron-lattice coupling giving rise to a quantitative description of the materials. In TTF-based systems, it is known that the amount of charge on the molecule can be probed by measuring the frequency of the central C=C bond vibration by Raman vibrational spectroscopy. Yue et al. (2010); Hassan et al. (2018) In fact, the low-temperature Raman spectrum of D-S shows new peaks at 1516 and 1405 cm*-1* that are not observed in the high-temperature phase. Ueda et al. (2014) These are assigned to the shorter and longer central C=C bonds of the w-D and w/o-D units, respectively. With a reliable description of the CO phase as demonstrated here, we expect that the phonon frequencies calculated on the basis of our results can now be directly compared with experiments. Tsumuraya et al. (2009)

The band gap is opened in our HSE06 calculations for D-S, and its value is in fairly good agreement with the experimental value estimated from the resistivity measurement. Ueda et al. (2014, 2015) On the other hand, for D-Se, our calculations exhibit a semi-metallic band structure while experiments show an insulating behavior. The full clarification of the accuracy of HSE06 for the evaluation of band gaps in strongly correlated molecular systems is left for future studies.

Finally let us discuss our results in light of previous theoretical studies for these compounds. Multi-component DFT based quantum chemistry calculations were performed by Yamamoto et al. Yamamoto et al. (2016, 2017) with a hybrid functional, the so-called Minnesota functional. The nuclear quantum effect of H/D on the potential energy surface of the H/D atoms was studied: It changes the potential surface from double well to single well. We consider that the two DFT-based calculation methods are complementary to each other: In our calculations, the quantum effect of H/D atom cannot be discussed, while the structural stability of the low-temperature phase in the actual crystalline form can be investigated in contrast to their calculations where only isolated molecular systems are treated.

On the other hand, Naka and Ishihara Naka and Ishihara (2018) studied an effective model where the correlated $\pi$-electron system is coupled with pseudospins describing the proton (H+) degree of freedom, assuming its double well potential. Their results show that two types of $\pi$-electron-proton coupled order compete with the dimer-Mott insulating state, both resulting from the cooperative effects of the proton ordering and CO. The two ordered states corresponds to the centrosymmetric and noncentrosymmetric phases discussed above. They considered the proton-electron coupling as an attractive potential, which is consistent with our results that the w-D monomers show lower HOMO energy than the w/o-D monomers (Fig. 3). However, in their model, the difference of inter-monomer distances in each dimer is not taken into account. Our results suggest that the electronic structure of the low-temperature phase is stabilized by both effects, namely, the variation of the intermolecular distances in two different dimers and the difference in MO energy levels depending on the position of the D atom. By the quantitative evaluation using first-principles calculations, these two effects indeed have the same order, as seen from our energy diagram in Fig. 3.

VI Summary

We studied the electronic and structural properties of the deuterium and $\pi$-electron coupled CO state of $\kappa$-D3(Cat-EDT-TTF)2 (D-S) and $\kappa$-D3(Cat-EDT-ST)2 (D-Se) using GGA and a range-separated hybrid functional HSE06. Using the experimental crystal structure, both GGA and HSE06 calculations show charge ordering but with a larger degree of charge disproportion in the latter, closer to the experimental situation. Using HSE06, an insulating band structure of D-S is obtained. By performing structural optimization with HSE06, the central C=C bond lengths that are sensitive to the degree of charge disproportionation are well reproduced, whereas GGA fails to stabilize the CO state. We also proposed possible patterns of D ordering and showed a stable noncentrosymmetric CO phase which has energy close to the experimentally realized phase.

VII ACKNOWLEDGEMENTS

The authors thank A. Ueda, H. Mori, M. Naka, H. Watanabe, T. Sasaki, and R. Kato for stimulating discussions. This research was funded by JSPS Grant-in-Aid for Scientific Research Nos. 16K17756, 26400377, 19K04988, 19K03723, 19K21860 and JST CREST Grant No. JPMJCR18I2. One the authors (T.T.) is partially supported by MEXT Japan, Leading Initiative for Excellent Young Researchers (LEADER). Computational work was performed under the Inter-university Cooperative Research Program and the Supercomputing Consortium for Computational Materials Science of the Institute for Materials Research (IMR), Tohoku University (Proposal No. 18K0090 and 19K0043). The computations were mainly carried out using the computer facilities of ITO at the Research Institute for Information Technology, Kyushu University, MASAMUNE at IMR, Tohoku University, NIMS Material Simulator, and HOKUSAI-GreatWave at RIKEN.