Comprehensive $ab$ $initio$ investigation of the phase diagram of quasi-one-dimensional molecular solids

Introduction—. Molecular solids serve as textbook materials to study many-body physics in condensed matter, despite their complex crystal structures with many atoms in the unit cell, owing to the success in modeling the low-energy electronic properties based on the frontier molecular orbitals [1, 2]. A general challenge is to understand multiple symmetry-broken phases observed in different systems in a systematic way, where the charge, spin, and lattice degrees of freedom are coupled to each other, based on the effective Hamiltonians and identifying the factors governing them [3].

The family of compounds constituted by TMTSF (= tetramethyltetraselenafulvalene) or TMTTF (= tetramethyltetrathiafulvalene) molecules, which we call together TM here, forming 2:1 salts with monovalent anions X^-1, TM₂X, has been studied for decades as a typical example showing such a rich variety of phases [4, 3, 2, 5]. TMTSF₂X historically outstands as the first organic system undergoing superconductivity [6]. Unconventional superconductivity with an anisotropic gap function is realized, whose ‘glue’ for the Cooper pair formation is the strong spin fluctuation, supported by the pressure-temperature phase diagram with the superconducting phase neighboring an antiferromagnetic (AFM) phase driven by the nesting of the Fermi surface [7, 8, 9, 10, 11].

Later, the isostructural TMTTF counterparts became realized as strongly-correlated electron systems, owing to the atomic substitution resulting in smaller bandwidths than in the TMTSF compounds [12, 5, 13]. In particular, at high temperatures of the order of 100 K, TMTTF₂X hosts a competition between a Mott insulating state, the so-called dimer-Mott (DM) insulator, and a charge ordering (CO) state with electron-rich and poor molecules spontaneously arranging to avoid the Coulomb repulsion: a common framework for the understanding of correlated quarter-filled systems [14, 15, 16]. In the low-temperature region ($\sim 10$K), on the other hand, where the magnetic interactions start to develop between the localized spins in the strongly correlated regime, the AFM state competes with the spin-Peierls instability accompanying lattice distortion owing to the quasi-one-dimensionality [17, 18, 19].

The variation among the symmetry-breaking states, in which both the charge and spin degrees of freedom are involved, depending on the applied pressure to the compounds and on the choice of the combination of TM and X (chemical pressure), has been extensively discussed [20, 21, 22, 23, 24]. Nevertheless, there is room toward a comprehensive understanding of them. There are different views on how the choice of X and the amount of applied pressure affect the microscopic parameters, and the impact of the difference between TMTSF and TMTTF is still obscure. Conventionally, the two systems are drawn in a unified pressure-temperature phase diagram [4, 5, 22, 23], as shown in Fig. 1(a). Although there are different versions in the literature, here we simply draw the ‘ambient pressure’ positions of the representative compounds in an equally spaced manner, together with the reported transition temperatures [27, 28, 29, 30, 13].

From the early days, in the TMTSF salts, it is discussed that the applied pressure increases the dimensionality [6]. On the other hand, the degree of dimerization along the one-dimensional chain direction was also pointed out to be an important factor, especially for the TMTTF salts [31]. And more recently, to understand the complicated interplay between the phases toward the left-hand end of the phase diagram, it was proposed that the CO instability can be the controlling factor [22]; this suggests the essential role of electron correlation. Whether we can choose a parameter that varies monotonically as we go along the compounds listed in the horizontal axis is the question we ask here.

In this work, to investigate these issues, we apply first-principles density functional theory (DFT) calculations [32, 33] to a number of members of the TM₂X family with centrosymmetric $X$; we numerically derive not only the transfer integrals between the molecules but also the correlation parameters. We then use these parameters as inputs to a highly-accurate numerical solver, i.e., the many-variable variational Monte Carlo (mVMC) method [34, 35, 36], for the effective extended Hubbard-type Hamiltonian and seek possible long-range ordering in the charge and spin degree of freedoms. As a result, we can directly compare with the experimentally observed phases and comprehensively extract the determining factor at a quantitative level.

Ab initio derivation of microscopic parameters—. The DFT calculations are performed using Quantum Espresso (version 6.6) [37] for the experimental crystal structures [13]. We employ norm-conserving pseudopotentials based on the Vanderbilt formalism with plane-wave basis sets [38, 39]. The exchange-correlation functional used in this study is the generalized gradient approximation by Perdew, Burke, and Ernzerhof [40], which is commonly used for accurately describing electronic states of molecular compounds [41]. The cutoff energies for plane waves and charge densities are $70$ and $280$ Ry, respectively. A $7\times 7\times 3$ uniform $\bm{k}$-point mesh was used with a Gaussian smearing method during self-consistent loops.

For deriving the interaction parameters, we use the constrained random phase approximation (cRPA) method [42, 43]. In fact, previous studies applying the cRPA method to molecular solids show good agreement between the results for the derived effective Hamiltonians and the experiments  [44, 45, 46, 47, 48, 49]. On the basis of the obtained DFT electronic states, we construct maximally localized Wannier function(MLWF)s and derive the parameters using RESPACK [50]. In the calculations, the energy cutoff for the dielectric function was set to be $3$ Ry.

Figures 1(b) and (c) show the band structure and the MLWFs, respectively, of TMTTF₂PF₆ as an example. The two bands crossing the Fermi energy are constituted from the bonding and anti-bonding combination of the highest occupied molecular orbitals situated on the two TMTTF molecules in the unit cell shown in Fig. 1(b). Using these MLWFs, we evaluated the transfer integrals and the density-density interactions, i.e. the on-site and off-site Coulomb interactions (the exchange interactions are negligibly small so we omit them in the following). The obtained microscopic parameters are listed in Table II in the supplemental materials [25]. We adopt the indices of bonds from refs. [22, 51], as shown in Fig. 1(d).

Then the following ab initio extended Hubbard-type Hamiltonian is obtained:

where $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ are the creation and annihilation operators of an electron with spin $\sigma$ at the $i$-th site, respectively. The number operators are defined as $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ and $N_{i}=n_{i\uparrow}+n_{i\downarrow}$.$t_{ij}$, $U$, and $V_{ij}$ represent the transfer integrals, on-site and off-site Coulomb interactions, respectively. We note that, although TM₂X show quasi-one-dimensional electronic band structures, the Coulomb parameters show large values also for the interchain bonds.

First, let us discuss the material dependence of intermolecular transfer energies based on our evaluation. Since it is known that the temperature variation of structural parameters is large [52], here we show the results for TMTTF₂X with X= NbF₆, AsF₆, Br, and PF₆ at fixed $200$ K, and for TMTSF₂X with X= NbF₆, AsF₆, and PF₆ at room temperature [25]. In Fig. 2, we show the variation of three different values along the horizontal axis of the phase diagram Fig. 1(a). The first two, the degree of interchain couplings and the degree of dimerization plotted in Figs. 2(a) and 2(b), respectively, have been discussed in the literature as mentioned above. We renormalize them by the mean transfer integral along the $a$-axis, $\bar{t}_{a}\equiv|(t_{a1}+t_{a2})/2|$, which is larger than those along the other directions.

The largest interchain transfer integral, which is $t_{b}$, is plotted in Fig. 2(a); the interlayer transfer integrals along the $c$-axis are negligible. Our results show larger $|t_{b}/\bar{t}_{a}|$ for the TMTSF salts than the TMTTF salts as a general trend; TMTSF has a broader extent of the wave function in the transverse direction than TMTTF, owing to the difference between the Se and S atoms at the sides of the molecules. As for the degree of dimerization, $\delta t_{a}/\bar{t}_{a}\equiv|t_{a2}-t_{a1}|/\bar{t}_{a}$ [Fig. 2(b)], the overall tendency is that it decreases from the left to the right side, while some deviates from the monotonic variation; for example, TMTTF₂Br shows a similar value to the TMTSF salts.

Finally, Fig. 2 (c) shows the on-site Coulomb interaction, $U$, which is a measure of the strength of electron correlation. The intermolecular parameters, $V_{ij}$, also show similar tendencies. We can see that $U/\bar{t}_{a}$ monotonically decrease from the left to the right side, except for the two salts at the right end: TMTSF₂PF₆ ($U/\bar{t}_{a}=9.59$) is slightly larger than TMTSF₂AsF₆ ($U/\bar{t}_{a}=9.46$). Since the electron correlation is responsible for all the phases, from the strongly correlated DM insulating and CO phases, to the nesting-driven AFM and the spin-fluctuation mediated superconducting phases, this behavior is consistent with the temperature scale for their occurrence, decreasing from left to right of the phase diagram in Fig. 1(a).

This tendency in $U/\bar{t}_{a}$ can be understood by the effect of chemical pressure. More concretely, first, by the unit cell volume decreasing from left to right, the transfer integrals reflecting the overlap between MLWFs increase. Next, since the screening effect becomes stronger with increasing the band width, the screened $U$ is smaller with increasing the pressure. As a combination of them, together with the difference in the bare onsite Coulomb parameter $U_{\rm bare}$ between the molecules, $U/\bar{t}_{a}$ monotonically decreases from left to right in Fig. 2 (c). For more detailed information about the parameters and their numerical accuracy, see ref. [25].

mVMC analysis—. To directly detect the symmetry-breaking phase transitions, we solve the $ab$ $initio$ effective Hamiltonian by using the mVMC method [34, 35, 36]. The form of the variational wave function used in this study is given by

where $\mathcal{P}_{\rm G}$ ($\mathcal{P}_{\rm J}$) represents the Gutzwiller factor [53] (Jastrow factor [54]) and $\mathcal{L}_{S}$ is the total spin projection. We use the spin singlet total spin projection ($S=0$) for the ground states. The pair-product wave function $|\phi_{\rm pair}\rangle$ is defined as

where $f_{ij}$ represents the variational parameters, and $N_{\text{s}}$ and $N_{e}$ are the number of sites and electrons, respectively. In the mVMC calculations, we map the lattice structure shown in Fig. 1(d) onto the square lattice without loss of generality. The total number of sites is given by $N_{\text{s}}=L_{x}\times L_{y}$. We impose a $4\times 2$ sublattice structure and the periodic boundary conditions, and optimize all the variational parameters simultaneously using the stochastic reconfiguration method [55, 34] to obtain the ground state.

We consider all the bonds that are listed in Table II in the supplemental materials [25]; they contain all the bonds shown in Fig. 1(d) as well as $V_{2a}$ that is the next-nearest-neighbor Coulomb interaction along the chain direction. We note that the effective Coulomb interactions are reduced by a constant shift $\Delta_{\rm DDF}=0.20$ eV [25], adopting the results of the previous study [56, 57], to take into account dimensional screening effects. We calculate the spin/charge structure factors defined as

where $N_{i}$ ($\bm{S}_{i}$) represents the number (spin) operator at site $i$. The Fourier transformation is performed for the mapped square lattice.

Let us present $S(\bm{q})$ and $N(\bm{q})$ for the two end members of the phase diagram, TMTTF₂NbF₆ and TMTSF₂PF₆. As shown in Figs. 3(a) and 3(b), for TMTTF₂NbF₆, which has the strongest electronic correlation, the spin (charge) structure factors has a sharp peak at $\bm{q}_{\rm peak}^{s}=(\pi/2,\pi)$ [$\bm{q}_{\rm peak}^{c}=(\pi,0)$]. The corresponding ordering pattern in the real space is shown in the inset of Fig. 3(a). Since these peak values remain finite after size extrapolation, the CO state with spin ordering is the ground state, consistent with previous studies [14, 22]. Note that, in the previous studies the long range Coulomb interactions are introduced by hand to stabilize the CO state with the ferroelectric pattern which is observed experimentally. We emphasize that, here the condition favoring the ferroelectric CO state is satisfied from the $ab$ $initio$ way, displayed in Fig. S1 (d) in the supplemental materials [25]. In contrast, we find that the peak values of the structure factors for parameters of TMTSF₂PF₆ become significantly small as shown in Figs. 3(c) and 3(d). This result is consistent with our results above that TMTSF₂PF₆ has much weaker electronic correlation.

To clarify the material dependence, we perform the mVMC calculations for all the members and conduct the size extrapolation for the obtained peak values, in the spin and charge structure factors, as shown in Figs. 4(a) and (b), respectively. The thermodynamic limits are defined as $M=\lim_{N_{\rm s}\rightarrow\infty}S(\bm{q}_{\rm peak}^{s})/N_{\rm s}$, $C=\lim_{N_{\rm s}\rightarrow\infty}N(\bm{q}_{\rm peak}^{c})/N_{\rm s}$. Both the spin and charge structure factors show nearly linear dependencies as a function of the linear dimension of the number of sites, i.e., $N_{\rm s}^{1/2}$.

The material dependence of $M$ and $C$ is shown in Figs. 4 (c) and (d), respectively, as a function of the inverse of the normalized Coulomb interaction, $(U-\Delta_{\rm DDF})/\bar{t}_{a}$. We find that the amplitude of the charge and spin orders decrease monotonically with respect to the amplitude of the Coulomb interactions. This tendency is consistent with the experimental results where the transition temperature of CO decreases by varying the material from left to right. At the border, TMTSF₂AsF₆ and TMTSF₂PF₆ show a small but finite value of $M$ of the order of 10^-2, and $C$ almost vanishes.

These results, combined with the ab initio results in Fig. 2(c), both clearly show that the charge and spin orderings are governed by the relative amplitude of the Coulomb interactions. This is the first $ab$ $initio$ evidence that the correlation effects play key roles in stabilizing the spin/charge ordering in TM${}_{2}X$ and serve as the good index for the horizontal axis of the unified phase diagram. The success here also indicates that our low-energy effective model is reliable at a quantitative level, considering the high precision of the mVMC method. This is owing to the fact that the band structure near the Fermi energy is simple, as shown in Fig. 1 (b), which guarantees the accuracy of the cRPA method, since it incorporates the effect of bands other than the target bands. In fact, in many of the molecular solids, bands near the Fermi energy are isolated, thus it is expected that our approach can be applicable to a wide range of functional molecular materials.

Perspectives—. We have shown that the material dependence in the unified phase diagram for TM${}_{2}X$ is governed by the correlation effect, represented by the on-site Coulomb interaction $U$, which gives a comprehensive understanding of the spin and charge orderings in this system. One point that disagrees with the experiments is that in our results, spin patterns with long periodicity such as incommensurate AFM (spin-density-wave) states suggested experimentally are difficult to reproduce in our finite-size systems. Finally, superconductivity is not found in our study, which is consistent with experiments showing it under pressure for the compounds we studied; its ab initio investigation is left for future studies.

Supplemental Materials for “Comprehensive $ab$ $initio$ investigation of the phase diagram of quasi-one-dimensional molecular solids”