Electronic correlation and geometrical frustration in molecular solids: A systematic ab initio study of

We systematically derive low-energy effective Hamiltonians for molecular solids β (cid:2) - X [Pd(dmit) 2 ] 2 ( X represents a cation) using ab initio density functional theory calculations and clarify how the cation controls the interdimer transfer integrals and the interaction parameters. The effective models are solved using the exact diagonalization method and the antiferromagnetic ordered moment is shown to be signiﬁcantly suppressed around the spin-liquid candidate of X = EtMe 3 Sb, which is reported in experiments. We also show that both the geometrical frustration and the off-site interactions play essential roles in the suppression of antiferromagnetic ordering. This systematic derivation and analysis of the low-energy effective Hamiltonians offers a ﬁrm basis to clarify the nature of the quantum spin liquid found in β (cid:2) -EtMe 3 Sb[Pd(dmit) 2 ] 2 . DOI

Introduction.Quantum spin liquids (QSLs) [1], which are Mott insulators without any broken symmetry, even at zero temperature, are new states of matter that have attracted much interest in the last few decades.In pioneering work by Anderson and Fazekas [2,3], it was proposed that geometrical frustration in the magnetic interactions can melt magnetic order and induce a QSL.Motivated by such a proposal, much theoretical and experimental work has been done to search for QSLs in frustrated magnetic materials [1,4,5].
Among several QSLs found in molecular solids, Pd(dmit) 2 salts offer an ideal platform for examining the key parameters that induce the QSL because it is possible to tune the ground states from the ordered states (antiferromagnetic ordering or charge ordering) to a QSL by systematically changing the cations [6].
In addition to that, because high-quality samples with less impurities are available, a detailed measurement of the thermal transport was done to clarify the nature of the QSL [12], and it was proposed that the large thermal conductivity indicates the emergence of an exotic particle such as the spinon in the QSL.However, this result has been challenged by recent experiments and the existence/absence of an exotic particle in the QSL is under hot debate [13][14][15][16].
Although theoretical studies such as those establishing lowenergy Hamiltonians are expected to play an important role for resolving the contradiction, there are few theoretical studies based on nonempirical methods.
In previous studies, a half-filled Hubbard model on an anisotropic triangular lattice was obtained as an effective microscopic Hamiltonian to describe the electronic structures in Pd(dmit) 2 salts [17][18][19] because the bands crossing the Fermi level are half filled and isolated from the other bands.These bands mainly originate from the antibonding pair of the highest occupied molecular orbitals (HOMOs) of the two Pd(dmit) 2 molecules that form a dimer [6,[20][21][22].A one-band model with three interdimer transfer integrals reproduces the density functional theory (DFT) bands well [17,23].Through extended Hückel calculations and tight-binding fitting to DFT bands, the amplitudes of the geometrical frustration, i.e., the anisotropy of the transfer integrals, have been evaluated.It has thus been proposed that the geometrical frustration governs the magnetic properties of the dmit salts and that the magnetic order can be suppressed by changing the transfer integrals [6].
In contrast to the transfer integrals, information on the interaction parameters is limited, even though they play an important role in the stabilization and suppression of magnetic order.Although the interaction parameters have been evaluated only for the QSL compound X = EtMe 3 Sb [17] in the studies of the dmit salts, the compound dependence and role in inducing the QSL state have yet to be clarified.We note that previous studies were limited to deriving low-energy effective Hamiltonians and it was not examined whether the low-energy effective Hamiltonians reproduce the cation dependence of the ground states in Pd(dmit) 2 salts by solving them.
In this Rapid Communication, to clarify the microscopic origin of the cation dependence of the ground states in Pd(dmit) 2 salts, we perform systematic ab initio derivations of the low-energy effective Hamiltonians, including both the transfer integrals and the interaction parameters for β -X [Pd(dmit) 2 ] 2 .At present, nine different Pd(dmit) 2 salts with a β -type structure are synthesized because the monovalent cations X can take three different types, i.e., MeY 4 , EtMe 2 Y 2 , and EtMe 3 Y , where the choice of pnictogen Y is P, As, and Sb.Through a comparison with the obtained lowenergy effective models, we find two trends in the parameters of the effective Hamiltonians: (1) Hopping parameters t c /t a [definitions are given in Fig. 1(a)] increase in the order of P, As, Sb. (2) On-site Coulomb interactions U/t a increase in the order of P, As, Sb (EtMe 3 Sb is exceptional).
Furthermore, to identify how the microscopic parameters affect the magnetic properties, the obtained models are solved using exact diagonalization [24], which enables us to clarify the overall trend of the magnetic properties of the frustrated magnets [25,26].As a result, the trend in the compound dependence of the magnetic ordered moment is successfully reproduced and the magnetic ordered moment is found to be significantly suppressed in X = EtMe 3 Sb, wherein the QSL state was reported experimentally.EtMe 3 Sb is sandwiched between the Néel-type antiferromagnetic order and the striped antiferromagnetic order when the possibility of charge ordering is ignored.We also show that both geometrical frustration and off-site interactions play key roles in the suppression of antiferromagnetic order.The present results clarify the microscopic origin of the QSL and offer a firm basis to comprehensively understand QSLs found in molecular solids.
Ab initio derivation of effective models.Based on ab initio calculations, the following single-band extended Hubbardtype Hamiltonian was obtained, where c † iσ (c iσ ) is a creation (annihilation) operator of an electron with spin σ in the Wannier orbital localized at the ith dmit dimers.The number operators are defined as n iσ = c † iσ c iσ and N i = n i↑ + n i↓ .We evaluate the transfer integrals t i j , the on-site Coulomb interaction U , the off-site Coulomb interaction V i j , and direct exchange interactions J i j in an ab initio way.We note that the double-counting problem on the Hartree terms for multiorbital systems [27,28] does not exist in this study because the employed model is a single-band model.From here, we detail how to evaluate these parameters.
To derive low-energy effective Hamiltonians, we first perform non-spin-polarized calculations with a first-principles DFT method [30,31] and obtain the Bloch functions for the experimental structures of the nine members of β -type X [Pd(dmit) 2 ] 2 , whose structures are measured above the antiferromagnetic or charge ordering transition temperature [32,33].All the structures have the same space group symmetry of C2/c.The positions of all the H atoms were relaxed because those determined from x-ray diffraction measurements typically show slightly shorter C-H bond distances than the DFT optimized positions.One-electron Kohn-Sham equations are solved self-consistently using a pseudopotential technique with plane-wave basis sets, which is implemented in the scalar relativistic code of QUANTUM ESPRESSO 6.3 [34].The exchange-correlation functional used is the generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof (PBE) [35].The cutoff energies for plane waves and charge densities were set to be 70 and 280 Ry, respectively.A 5 × 5 × 3 uniform k-point mesh was used with a Gaussian smearing method during the self-consistent loops.
After obtaining the Bloch functions, the maximally localized Wannier functions (MLWFs) were constructed using RESPACK [36].To make the MLWFs, the half-filled bands crossing the Fermi level were selected as the low-energy degrees of freedom.The initial coordinates of the MLWFs were set to be at the center between two [Pd(dmit) 2 ] monomers to generate a one-band model (the so-called dimer model) [37].
In Fig. 1(a), we show the MLWF for X = EtMe 3 Sb.The MLWF, whose center position is located at the center between two [Pd(dmit) 2 ] monomers, spreads over the molecule and forms the dimer unit.Figure 1(b) shows that the MLWF nearly perfectly reproduces the band structures obtained by the DFT calculations.Using the MLWF, the transfer integrals of the low-energy effective models were evaluated as where φ n,R is the nth MLWF centered at R, and H k is the onebody part of the ab initio Hamiltonian.The lattice structure in the dimer units is shown in the inset of Fig. 1  that t a is the largest transfer integral and t c /t a denotes the amplitude of the geometrical frustration as discussed later.
Figure 2(a) shows the compound dependence of the normalized transfer integrals obtained by the MLWF fitting.For comparison, the transfer integrals obtained by extended Hückel calculations are shown.Although the transfer integrals evaluated by the MLWF are slightly different from those obtained by extended Hückel calculations, the trend of the compound dependence is consistent, i.e., t b /t a is not largely dependent on the cations and t c /t a increases in the order of P, As, and Sb.This trend is also found in the literature and the origin can be attributed to the cation radius, which controls the distortion of the dmit molecules [32].
The interactions were also evaluated by the constrained random-phase approximation (cRPA) [38] method using RESPACK.The energy cutoff for the dielectric function was set to be 3 Ry.The interaction terms are given as follows, where H W represents the interaction term of the ab initio Hamiltonians.Although three-body and four-body interactions generally occur, we only treat the two-body interactions, such as density-density interactions U mn (R) = W mm,nn (0, 0, R, R) (that is, the on-site and off-site Coulomb interactions) and the direct exchange interactions J mn (R) = W mn,nm (0, R, R, 0) because the amplitudes of the other terms are negligibly small.In the Supplemental Material (SM) [39], we show the dependence of the number of screening bands on the interaction parameters.We note that when the number of screening bands becomes large, the cRPA method can be valid in contrast to previous studies [40,41].Figure 2(b) shows the compound dependence of the normalized on-site Coulomb interaction U/t a .U/t a increases in the order of P, As, and Sb, as in the transfer integrals, except for EtMe 3 Sb.As we have detailed in the SM [39], the compound dependence of U/t a is mainly controlled by changes in t a because U is not largely dependent on the compounds.The transfer integral t a decreases in the order of P, As, and Sb because the distortions of the molecules become larger in this order [32].In contrast, U does not show a large compound dependence (see S. 1 in the SM [39]), and for this reason, U/t a increases in the order of P, As, and Sb.Low-temperature structures were employed for EtMe 3 Sb; therefore, t a becomes large and U/t a becomes small.This is the origin of the exceptional behavior observed in EtMe 3 Sb.We also point out that U/W (W is the bandwidth) is estimated as U/W ∼ 2.0 for EtMe 3 Sb, which is roughly consistent with the experimental estimation (U/W ∼ 2.3) [42].
Analysis of effective models.To clarify how the geometrical frustration and the interaction parameters affect the magnetic properties in the low-energy effective models defined in Eq. ( 1), exact diagonalization was performed for small clusters (the system size is N s = 4 × 4 with periodic boundary conditions).In the calculations, we take the transfer integrals, off-site Coulomb interactions, and the direct exchange interactions up to the next-nearest neighbor.Using exact diagonalization, the compound dependence of the magnetic properties can be clarified without relying on any specific approximations.It should be noted that it is necessary to introduce a constant shift in the interaction parameters to reflect the two-dimensionality of the effective models, and to obtain physically reasonable results.Following the parameter-free approach, we employ a constant shift of = 0.30 eV for all compounds, which is a value comparable to that in the literature [17].In S. 2 in the SM [39], we examine the effects of and confirm that the changes in the constant shift do not change the results significantly.
Figure 3(a) shows the spin structure factors, for several compounds, i.e., X = Me 4 P, EtMe 3 Sb, and Et 2 Me 2 Sb (the spin structure factors for all nine compounds are shown in S. 3 in the SM [39]).In Me 4 P, which shows the highest Néel temperature and anisotropy of transfer integrals, the spin structure factor is a sharp peak at q = (π, 0), which indicates a stripe magnetic order [top panel in Fig. 3(b)].
The amplitudes of the spin structure factors are significantly suppressed in EtMe 3 Sb and no clear signature of magnetic order is evident.Additionally, we note that no clear signatures of other exotic nonmagnetic phases, such as the bond-order phase [43,44], are observed (see S. 4 in the SM [39]).This indicates that the QSL state is formed in this compound.On the other hand, for Et 2 Me 2 Sb, a sharp peak appears at q = (π, π ), which indicates a Néel-type magnetic order shown in the bottom panel in Fig. 3(b).The ground state is a charge ordered state [45]; therefore, this result is apparently inconsistent with previous experimental results.However, this discrepancy can be attributed to the discarding of the longrange part of the Coulomb interactions; in a 4 × 4 system size, only the Coulomb interactions up to the next-nearest neighbor can be treated, and the charge order pattern cannot be accurately determined.We note that the interplay of longrange Coulomb interactions and electron-phonon couplings stabilizes the charge ordered phase in Et 2 Me 2 Sb [46].Our analysis indicates that Et 2 Me 2 Sb has an antiferromagnetic instability toward (π, π ), if we ignore the instability toward charge ordering.
Figure 3(c) shows the compound dependence of the peak values of the spin structure factor S(q peak ) as a function of t c /t a .X = EtMe 3 Sb is located around the boundary of the stripe and Néel magnetic order, and S(q peak ) is significantly reduced at X = EtMe 3 Sb.The overall trend of the compound dependence of the spin structure factors is consistent with experimental results [6].It is especially noteworthy that the low-energy effective Hamiltonian for EtMe 3 Sb shows a suppression of the spin structure factors, which is consistent with the experimentally observed QSL behavior.
To clarify the origin of the magnetic ordered moment reduction, we systematically changed the parameters in the  ab initio effective Hamiltonian for EtMe 3 Sb.First, the effects of the geometrical frustration were examined by changing t c /t a .In Fig. 4(a), we show the t c /t a dependence of S(q peak ).By artificially decreasing t c /t a , a sudden change in S(q peak ) occurs at t c /t a ∼ 0.58, which indicates a first-order phase transition between the stripe magnetic ordered phase and the possible QSL state.This result clearly shows that the geometrical frustration plays a key role in the suppression of the magnetic ordered moment.

S(q)/Ns
Next, we introduce λ, which monotonically scales the offsite interactions including the off-site Coulomb interactions (V i j ) and the direct exchange interactions (J i j ), i.e., (V i j , J i j ) is scaled as (λV i j , λJ i j ).Note that λ = 1 corresponds to the ab initio Hamiltonian and λ = 0 corresponds to the simple Hubbard model that has only on-site Coulomb interactions U .The dependence of S(q peak ) on λ is shown in Fig. 4(b).S(q peak ) increases by decreasing λ, and the stripe magnetic order appeared below λ = 0.5 (the spin structure factor at λ = 0 is shown in the inset).This result indicates that the offsite interactions, which were often ignored in previous studies, suppress the magnetic ordered moment and play an important role in the stabilization of the QSL state in Pd(dmit) 2 salts.We note that both V i j and J i j reduce the magnetic ordered moment (see S. 5 in the SM [39]).The mechanism of the reduction can be attributed to the reduction of U and resultant enhancement of the higher-order interactions such as the ring-exchange interactions, which suppress the magnetic long-range orders [47][48][49].
Summary.To conclude, we have successfully reproduced the overall trend in the magnetic properties by performing comprehensive ab initio calculations for all available β -type Pd(dmit) 2 salts.An antiferromagnetic order with a stripe pattern becomes the ground state for small t c /t a such as Me 4 P.This result is consistent with experimental results that the Néel temperatures increase with decreasing t c /t a .However, the magnetic ordered moment is significantly suppressed around EtMe 3 Sb.These results indicate that a QSL state without magnetic order appears in EtMe 3 Sb, which is consistent with previous experimental results.We have also shown that cooperation of the geometrical frustration and the off-site interactions induces the suppression of the magnetic ordered moment in EtMe 3 Sb.To analyze the QSL state in molecular solids, most theoretical calculations have been conducted for frustrated Hubbard models that only have on-site Coulomb interactions [19,50,51].The present ab initio calculations, which suggest the importance of off-site interactions, require reconsideration of the use of such simple models to describe the QSL state in Pd(dmit) 2 salts.
The analysis of the low-energy effective model in this Rapid Communication is limited to small system sizes in order to perform an exact analysis to clarify the general trend in the Pd(dmit) 2 salts.Several highly accurate wave-function methods for strongly correlated electron systems have recently been developed [52][53][54].A clarification of the possible exotic elementary excitation of the QSL state in EtMe 3 Sb using such cutting edge methods, which is a hotly debated issue in experiments [12][13][14]16], is a significant challenge but is left for future studies.We also note that applications of the employed ab initio method to other families of molecular solids such as κ-(BEDT-TTF) 2 X salts will help us to comprehensively understand the nature of the QSL found in other compounds [55].

FIG. 1 .
FIG. 1.(a) Wannier function for EtMe 3 Sb drawn using VESTA [29].The inset shows a schematic of the lattice structure of EtMe 3 Sb in the a-b plane.Circles correspond to Wannier centers and transfer integrals between the nearest-neighbor Wannier orbitals are shown.(b) Band dispersion for EtMe 3 Sb.The dotted (solid) lines are obtained by the first-principles DFT method (Wannier interpolation).

FIG. 2 .
FIG. 2. (a) Band dispersion for EtMe 3 Sb.The dotted (solid) lines are obtained by the first-principles DFT method (Wannier interpolation).(b) Material dependency of t b /t a and t c /t a .The dotted (solid) lines are obtained by the transfer integrals calculated by the extended Hückel method [32] (Wannier function basis).

FIG. 3 .
FIG. 3. (a) Spin structure factors for three typical compounds.A significant reduction of the Bragg peak occurs in EtMe 3 Sb.(b) Schematic diagrams for stripe and Néel ordered states.(c) Dependence of the spin structure factors on t c /t a .For each family, t c /t a increases in the order of P, As, and Sb.

FIG. 4 .
FIG. 4. (a) Dependence of S(q peak ) on t c /t a .Around t c /t a ∼ 0.57, the stripe magnetic correlations suddenly increase.The inset shows the spin structure factors at t c /t a = 0.5, which indicate a stripe magnetic order.(b) Dependence of the peak value of the spin structure factors on λ for X = EtMe 3 Sb.The inset shows the spin structure factors at λ = 0.