Unconventional dual 1D–2D quantum spin liquid revealed by ab initio studies on organic solids family

Organic solids host various electronic phases. Especially, a milestone compound of organic solid, β′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta ^{\prime}$$\end{document}-X[Pd(dmit)2]2 with X=EtMe3Sb shows quantum spin-liquid (QSL) properties suggesting a novel state of matter. However, the nature of the QSL has been largely unknown. Here, we computationally study five compounds comprehensively with different X using 2D ab initio Hamiltonians and correctly reproduce the experimental phase diagram with antiferromagnetic order for X=Me4P, Me4As, Me4Sb, Et2Me2As and a QSL for X=EtMe3Sb without adjustable parameters. We find that the QSL for X=EtMe3Sb exhibits 1D nature characterized by algebraic decay of spin correlation along one direction, while exponential decay in the other direction, indicating dimensional reduction from 2D to 1D. The 1D nature indeed accounts for the experimental specific heat, thermal conductivity and magnetic susceptibility. The identified QSL, however, preserves 2D nature as well consistently with spin fractionalization into spinon with Dirac-like gapless excitations and reveals duality bridging the 1D and 2D QSLs.


INTRODUCTION
Organic solids offer plenty of playgrounds of conspicuous phenomena including superconductivity competing with charge and spin orders, and semi-metallic excitations with Dirac dispersions 1,2 . However, the severe competitions of the diverse phases and mechanisms of the phenomena found in complex crystal structures of the organic solids remain challenges of condensed matter physics.
Despite a large number of atoms in the unit cell of such organic compounds, however, the band structure near the Fermi level is mostly simple consisting of only LUMO (lowest unoccupied molecular orbital) and HOMO (highest occupied molecular orbital) around the Fermi level. Large inter-molecular distance leads to narrow bandwidths while poor screening of Coulomb interaction due to the sparse bands near the Fermi level leads to strongly correlated electron systems with various types of Mott insulators.
Among them, the QSL phase was proposed as novel states of matter in two families of correlation driven Mott insulating compounds, namely, κ-(ET) 2 X (ET= bis(ethylenedithio) tetrathiafulvalene) with the anion X=Cu 2 (CN) 3 3 and β 0 -X[Pd(dmit) 2 ] 2 (dmit=1,3-dithiole-2-thione-4,5-dithiolate) with the cation X=EtMe 3 Sb 4 , where apparent spontaneous symmetry breaking including the antiferromagnetic (AF) order is unusually absent even at several orders of magnitude lower temperatures than the spin-exchange interaction. Although the geometrical frustration arising from the triangular lattice structure of the dimerized Pd(dmit) 2 or ET molecule hosting a spin-1/2 electronic spin looks important for the absence of the AF order, its nature and mechanism of the emergence endorsed by a quantitative estimate are missing. In addition to the QSL, if one replaces X with other ions, it shows a variety of phases including the AF order, charge order and valence bond solid 1,2,5 . Because of severe competitions of these phases, ab initio approaches without adjustable parameters are desired particularly for the enigmatic QSL to reach realistic understanding.
We focus on a typical candidate of QSL, X[Pd(dmit) 2 ] 2 . In fact, it shows 4,6 smaller broadening of the nuclear magnetic resonance (NMR) spectra than κ-(ET) 2 Cu 2 (CN) 3 1,7 implying smaller effects of extrinsic spatial inhomogeneity. At low temperatures of the QSL material, EtMe 3 Sb[Pd(dmit) 2 ] 2 , the specific heat 8 and the thermal conductivity κ 9 are reported to be proportional to temperature T, although controversies exist for the latter [10][11][12] . The magnetic susceptibility stays at a nonzero constant 13 . These are roughly similar to the conventional Fermi-liquid metal but of course the electrons are frozen in the Mott insulator and the spin degrees of freedom must be responsible for them. In addition, the relaxation rate of NMR, 1/T 1 seems to be scaled by T 2 in the range 0.05K ≤ T < 1 K 6 in contrast to the conventional Fermi-liquid behavior ∝ T. It is then desired to gain insights into the experimental consistency from ab initio calculations.
Here, we thoroughly study purely ab initio electronic Hamiltonians without any adjustable parameters for five dmit compounds using the experimental structure without optimizing lattice structures except for positions of hydrogen atoms. This ab initio study allows us to correctly reproduce the overall experimental phase diagram at low temperatures of β 0 -X[Pd(dmit) 2 ] 2 with X=Me 4 P, Me 4 As, Me 4 Sb, Et 2 Me 2 As exhibiting the AF order and with X=EtMe 3 Sb exhibiting QSL phases. We do not consider β 0 -Et 2 Me 2 Sb[Pd(dmit) 2 ] 2 , because the observed charge order with the simultaneous spontaneous lattice distortion in this compound 14,15 requires structural optimization, which is out of the scope of this paper. Then we conclude that the QSL phase has a character of primarily 1D spin liquid but combined with 2D nature, which accounts for the experimental properties consistently. From the study, we attempt to pin down the nature of QSL observed in the experiments, which is one of the central issues in condensed matter physics. 1

Ab initio framework
Before going into our results, we summarize our framework of the calculation for the sake of readers to understand definitions of quantities, which we will analyze. We use ab initio low-energy effective Hamiltonian consisting of the half-filled HOMO band crossing the Fermi level derived for β 0 -X[Pd(dmit) 2 ] 2 [16][17][18] . For the procedure and the obtained parameters for the Hamiltonians, see Methods. It is built on the weakly coupled two-dimensional layers and the Pd(dmit) 2 dimer forms an anisotropic triangular lattice within a layer. We assign theã,b andc bonds, for which the Hamiltonian presented below (Eq. (1)) contains the transfer parameters t a , t b and t c , respectively, as is illustrated in the middle panel of Fig. 1. The bonds are chosen so that t a , t b and t c are ordered from the stronger to weaker amplitudes. In the actual calculation, we take a deformed lattice structure illustrated in the right panel of Fig. 1 just for computational simplicity. We sometimes use x and y axes for the plot of the 2D Brillouin zone and note that x and y correspond toã andc directions, respectively. (Do not confuse the crystallographic a, b and c axes with the present assignment ofã,b andc axes).
Because of very weak interlayer hopping (the largest interlayer hopping is less than 0.8 meV), the 2D plane is sufficient to capture physics of our interest about the ground state here. Then our Hamiltonian is built on the two-dimensional plane of the dmit salts. We take into account the effect of interlayer interaction through the dimensional downfolding method established in organic solids and iron-based superconductors 16 (see Methods). The maximally localized Wannier function of the HOMO orbital, whose band crosses the Fermi level, is constructed for the molecular orbital 19,20 . The lattice structure in the simulation is depicted later in Fig. 6, which is, despite the deformation of the shape, able to sufficiently describe the ab initio Hamiltonian as the network of the Pd(dmit) 2 dimers, where a Wannier orbital is extended over a dimer.
The resultant ab initio single-band effective Hamiltonian derived in refs. 17,18 has the form where i, j represent the dimer indices, and c y iσ (c iσ ) is the creation (annihilation) operator of electrons with spin σ (=↑ or ↓) at the ith Wannier orbital, and the number operator is n i = ∑ σ n iσ with n iσ ¼ c y iσ c iσ . Here, t ij is the hopping parameters depending on the relative coordinate vector r i −r j , where r i is the position vector of the center of the ith Wannier orbital. In the present study, for t ij , we retain the nearest neighbor pair of r i and r j in each direction of a,b andc as is illustrated in Fig. 1, while V ij is retained up to the third neighbor (see Fig. 6). Here, U and V ij are the screened on-site and off-site Coulomb interactions, respectively, as is illustrated in Fig. 6. Note that the Hamiltonian parameters of both transfer and interaction terms contain neither adjustable parameters nor fitting and are determined solely by using the experimental lattice structure at low temperatures and by following the established procedure of the maximally localized Wannier functions. The spatial ranges of the transfer and the interaction, which we take are sufficient to reach the convergence for the present study and the values at longer distance are small. We also ignore the direct exchange interactions because they are at most 3 meV and we expect it does not alter our conclusions. See Methods for the parameter values of Hamiltonian (1) used in the present study.
We apply a variational Monte Carlo (VMC) method [21][22][23] to the ab initio Hamiltonians to reach highly accurate ground states. For details of the numerical method, see Methods. Figure 2 is one of our central results of this paper, showing the agreement between the experimental (bottom plane) and the calculated (vertical plane) material dependences of the lowtemperature phases. In the calculated results, the collinear AF and the QSL states are the two lowest energy states among severely competing various candidates all through the five compounds studied. Therefore, we plot the calculated energy difference between these two in the vertical plane. We find that the experimental QSL is reproduced in the calculated ground state for X=EtMe 3 Sb. Aside from a severe competition for X=Et 2 Me 2 As, which is indeed suggested by the experimental coexistence of AF and QSL, Fig. 2 shows that our ab initio results successfully reproduce the AF ground state observed in the experiments of X=Me 4 P, Me 4 As, and Me 4 Sb, as shown in the bottom plane. The AF state has the Bragg peak at (π, 0) in the notation of Fig. 1. The AF stability was recently studied by the quasi-1D approach based on the random phase approximation 24 .

Phase diagram
Note that the abscissa of Fig. 2 shows only (t c − t b )/t a dependence among the ab initio parameters, while the plots of the five materials are the results of computation by using the full ab initio parameters. The overall monotonic dependence of Δe shows that (t c − t b )/t a is indeed a principally important parameter for the evaluation of the phase diagram. The experimental Néel temperature in the bottom plane is also consistently ordered with (t c − t b )/t a and further supports the relevance and accuracy of the ab initio effective Hamiltonian, because it is natural to expect that the Néel temperature is linked to the relative stability of the ground state. Within the parameter control of (t c − t b )/t a , the transition between the AF and QSL phases is a clear first-order transition, represented by the level crossing. The importance of the full anisotropy of the transfer in the inequilateral triangular lattice to stabilize the QSL phase was pointed out by the where a Pd(dmit) 2 molecule is depicted as a long gray oval. b Modeled triangular lattice with three different electronic transfers. The strongest, middle and weakest transfers, t a (on the red bondã), t b (on the blue bondb), and t c (on the green bondc), respectively, (they correspond to t B , t S, and t r , respectively in the literature 54 ). c Deformed structure on a L × L lattice with the system size N s = L 2 used in the present calculation.
projected BCS study of the Hubbard-type model by using the ab initio hopping 25 .
Nature of quantum spin liquid and antiferromagnetic state Figure 3 shows the spin structure factor S(q) (the spin correlation in the Fourier space) for X=Me 4 P as a representative case of the AF ordered compounds and for the QSL compound, X=EtMe 3 Sb (see Methods for details of the calculation method and definition of physical quantities). For the AF state of X=Me 4 P, a strong sharp peak at q = (π, 0) indeed indicates the conventional collinear 2D AF order. The 2D ordering pattern is illustrated in Fig. 2. Other compounds, X=Me 4 As, Me 4 Sb, and Et 2 Me 2 As also show the same type of the order. In the QSL state for EtMe 3 Sb, S(q) also has a peak at the same momentum q = (π, 0), but with a substantially reduced height. Interestingly, a ridge line along q x = π emerges in S(q) of the QSL state, implying anisotropy between the x and y directions (namely between the chain (ã) and interchain (b) directions). Such anisotropy is in sharp contrast with the isotropic order in the AF phase. Furthermore, the feature of the anisotropy is largely different from the previous numerical studies for the anisotropic triangular Hubbard or Heisenberg model 26 . Figure 3(c) shows the size dependence of the peak value of S(q) divided by the system size N s = L 2 . Because the AF long-range order is signaled by a nonzero S(q)/N s in the thermodynamic limit 1/L → 0, AF order exists for X=Me 4 P and it is absent for X=EtMe 3 Sb. The spin-spin correlations in real space C(r) for the QSL state of X=EtMe 3 Sb shown in Fig. 3d with the log-log plot in the x (namely, t a ) direction indicates a power law C(r) ∝ |r| −p with p = 1.88 ± 0.01 suggesting the algebraic QSL with the gapless excitation. On the other hand, Fig. 3(e) shows a semilog plot of the correlations in the y direction, for which the exponential form CðrÞ / exp½Àjrj=ξ ? with the correlation length ξ ⊥~1 .33 ± 0.04 is appropriate and thus the excitation is gapped in the interchain (y) direction at least within the system size studied here. Of course, the same exponential decay of the correlation is observed in the t b (namelyb) and t c (namely,c) directions. The anisotropic correlation develops the ridge in the q y dependence of S(q).
Finally, the spin Drude weight D s calculated in the QSL state for X=EtMe 3 Sb is shown in Fig. 4. The size dependence of D s for the QSL shown in Fig. 4(b) assures that D s is nonzero and large in the thermodynamic limit, consistently with the power-law decay of the spin correlation in the x direction. In contrast, the spin Drude weight for the interchain y component is vanishing.
All the results support that the QSL with one-dimensional and gapless excitation is realized for X=EtMe 3 Sb, which is our second central result of the paper.

DISCUSSION
The QSL of β 0 -EtMe 3 Sb[Pd(dmit) 2 ] 2 is characterized by the gapless spin excitation and associated algebraic power-law decay of the spin correlation in the direction of the largest transfer t a (namely, along the chain direction). The power p~1.9 for CðrÞ / exp½iQ Á rr Àp is similar to the value obtained for the QSL in the squarelattice J 1 -J 2 Heisenberg model with the nearest (J 1 ) and the nextnearest-neighbor (J 2 ) exchange interactions, where p~1.4-1.7 27 . However, in contrast to isotropic 2D spin correlation in the J 1 -J 2 model, the correlation decays exponentially in the interchain direction with the correlation length ξ ⊥~1 .3 lattice constant. The anisotropy makes the peak height of S(q) non-divergent in the thermodynamic limit in contrast with the 2D J 1 -J 2 model.
One might regard the present QSL as essentially the same as the 1D spin liquid 26,28,29 smoothly connected to an effective 1D Heisenberg or Hubbard model. However, the revealed properties are not so simple. First, the nonzero correlation length in the interchain direction makes a prominent peak at (π, 0) in S(q) commonly to the case of the 2D J 1 -J 2 Heisenberg model. In the decoupled arrays of 1D Heisenberg chains, one would expect an equal-height ridge along the line (π, 0)-(π, π) with logarithmically divergent height where p = 1 and without q y dependence. In the present QSL, the ridge exists but the height is not divergent even at (π, 0) because p > 1. For the moment, we leave the issue about the possible presence or absence of the transition between the pure isolated chain and the present case for the future study. Substantial reduction of S(q) toward (π, π) from (π, 0), namely partially 2D-like correlation implies a small but nonzero dispersion of the excitation in the interchain direction, which may make an energetic hierarchy.
For the 2D J 1 -J 2 model, it was suggested that the gapless excitation is well accounted for by the fractionalization of a spin into two spin-1/2 spinons, where the spinon has the Dirac-like linear dispersion at (±π/2, ±π/2) 27,30 . Let us discuss the possibility of the fractionalization for the present QSL of β 0 -EtMe 3 Sb[Pd (dmit) 2 ] 2 . To gain insight into the possible existence of spinon and to elucidate its dispersion in the present QSL, we have analyzed the structure of our wavefunction for the QSL ground state of β 0 -EtMe 3 Sb[Pd(dmit) 2 ] 2 as is described in Methods. The result shows that the elementary excitation is consistently described by the spinon born out as a consequence of the fractionalization of a spin, where the spinon has Dirac-type gapless dispersion around (±π/2, 0) and (±π/2, π) as is schematically illustrated in Fig. 5. Because the measurable spin is constructed from the two-spinon excitation, the gapless points for the triplet excitations appear at around (0, 0), (±π, 0), (0, ±π), and (±π, ±π). This structure has an essential similarity to the 2D J 1 -J 2 Heisenberg model 27 , implying the 2D nature of the QSL.
However, the exponential decay of the spin correlation and vanishing spin Drude weight in the interchain direction suggest a very anisotropic and small or incoherent dispersion away from the

Ab initio Exp
nodal point in the interchain direction in contrast to the isotropic conic dispersion for the 2D J 1 -J 2 model, generating the QSL of roughly 1D-like nature. By considering the one-dimensional anisotropy, it is insightful to compare the present QSL with that of the 1D antiferromagnetic Heisenberg model defined by with J > 0 for the spin-half operator S i ¼ ðS x i ; S y i ; S z i Þ at site i. The specific heat of the 1D Heisenberg model is given by the T-linear 31 , where k B is Boltzmann constant and N A is Avogadro constant. If we employ the strong coupling expansion for the ab initio Hamiltonian, the superexchange interaction J a for the two neighboring spins in theã (namely, t a ) direction is estimated as J a ¼ 4t 2 a =ðU À V 1 Þ $ 0:031 eV for β 0 −EtMe 3 Sb[Pd(dmit) 2 ] 2 , by using the ab initio parameters listed in Table 1 of Methods. Then γ is estimated as γ~15 mJ K −2 The uniform magnetic susceptibility of the 1D Heisenberg model is B =ðJ a π 2 Þ with the g-factor and the Bohr magneton μ B 32,33 . By using J a~0 .031eV, and g~2, we obtain χ=μ 2 B $ 13 eV −1 , while the experimental value is χ=μ 2 B $ 12À24 eV −1 per the dimer 4,13 , which is consistent with each other by considering the experimental uncertainty.
Frequency dependent thermal conductivity κ of the 1D Heisenberg model is estimated to be κ 1D (ω) = κ s δ(ω) with κ s /T = π 3 J/(6aℏ 2 ), where a is the lattice constant 34 . If the delta function is replaced by the Lorentzian ω W =πðω 2 þ ω 2 W Þ to take into account the lifetime τ = 2π/ω W , we obtain κ 1D s ðω ¼ 0Þ=T ¼ π 2 J a =ð6a_ 2 ω W Þ. By considering a~10 −9 m and J a~0 .031eV, κ 1D s ðω ¼ 0Þ=T is estimated as 7.7 × 10 22  , the unit of the momentum q x and q y is π. c Size dependences of the peak value of S(q) denoted as S(q peak ) in the AF state for X=Me 4 P and the QSL state for X=EtMe 3 Sb. Dashed lines are guides for the eye. The AF state scales to a nonzero S(q peak )/N s in the thermodynamic limit, evidencing for the AF long-range order. On the other hand, extrapolation to zero in the thermodynamic limit for the QSL state confirms the nonmagnetic state. d Log-log plot of C(r = (r, 0)) for X=EtMe 3 Sb. Black thin line is the fitting curve C fit (r) ∝ ∑ n |nL + r| −p with n being integer by taking account of the antiperiodic boundary condition, which is fit by p = 1.88 for L = 24. We use data points for 2 < r < L − 2 when we obtain the fitting curve. e Semilog plot of C(r = (0, r)) for X=EtMe 3 Sb. Black thin line is the fitting curve C fit ðrÞ / P n expðÀjnL þ rj=ξ ? Þ ½ with ξ ⊥ = 1.33 for L = 24. Error bars are obtained from the Monte Carlo sampling and the variance extrapolation, details of which are described in Methods. The length unit is the lattice constant in panels (d) and (e). In panels (c-e), error bars obtained in the VMC calculation are plotted and if not they are smaller than the symbol size. the carrier relaxation time τ = 2π/ω W~9 .0 × 10 −12 s. A simple expectation of the spin velocity v s = J a a/ℏ results in the mean free path ℓ s = v s τ = 0.42 μm, which is a value comparable to the reported estimate~1 μm 9 . More precisely, the experimental value interpreted by the 1D Heisenberg model suggests ℏω W /J~0.015 implying small damping in the order of 10 −2 for the propagation induced plausibly by the spinon dynamics driven by the energy scale of J. Therefore, although controversies exist as is mentioned above, if a large thermal conductivity experimentally reported is intrinsic, it can be essentially accounted for by the interpretation of 1D-like QSL. The obtained ℓ s is substantially longer than that of the inorganic compound Cs 2 CuCl 4 35 , where a 1D QSL was claimed. The present result implies that the 1D QSL found here has potentially much longer τ. Then all of the above thermodynamic and transport properties can be interpreted by the essentially 1D-like QSL found here.
The spin-lattice relaxation time T 1 reported as scaled by 1/T 1 ∝ T 2 below 1K 6 , implying the point-like gapless triplet excitation looks contradicting the T-linear specific heat and nonzero magnetic susceptibility compatible with the constant density of states. However, a consistent picture may emerge, if a spinon with a highly anisotropic nodal and gapless dispersion exists at~(±π/2, 0) and~(±π/2, π), where the 2D nature could be detectable only at low temperatures below 1K and by measuring the gapless momenta with an appropriate form factor. The power of the dispersion around the node in the chain direction is an intriguing issue but is beyond the scope of the present paper.
It turned out that the Dirac-type gapless excitation for the 2D QSL can be connected to the vanishing spinon dispersion along the Fermi line for T > 1 K, which behaves as if it is essentially the 1D QSL. β 0 -EtMe 3 Sb[Pd(dmit) 2 ] 2 is located in this parameter space close to the 1D QSL. Furthermore, the controversy and sensitivity of the thermal conductivity [10][11][12] are well understood from the one-dimensional nature: The exponential decay of the spin correlation and the vanishing Drude weight imply the vanishing thermal conductivity in the interchain direction. Then even small angle misalignment of the single crystal either in samples or in measurements causes serious reduction of κ and serious sample or measurement dependence. Effects of randomness may also seriously disturb the ideal behavior of the 1D-like spin liquid 11 . Observed very sensitive dependence might be able to be interpreted from the 1D nature of the QSL, which has also been pointed out in ref. 36 (See also Methods).
In summary, we have studied the family of dmit organic salts by using the ab initio Hamiltonians of 5 compounds. The obtained low-temperature phases are consistent with the experimental reports as the AF state for X=Me 4 P, Me 4 As, Me 4 Sb, Et 2 Me 2 As and the QSL for EtMe 3 Sb. The relative stabilities of the four AF  Two red spinons, one at (−π/2, π) and the other at (−π/2, 0) in (a) can make measurable gapless triplet excitation shown in (b) at (−π, π). Two blue spinons at (π/2, 0) in (a) can make measurable gapless triplet excitation shown in (b) at (π, 0).
compounds increasing with decreasing t c − t b are correctly reproduced in the ab initio calculation in agreement with the experimental trend. Thanks to this firm correspondence, the nature of the QSL in the real material is identified as the 1D-like spin liquid: The spin correlation decays exponentially and the spin Drude weight vanishes in the interchain direction. Though a controversy exists in the experimental reports of the thermal conductivity, the specific heat, the susceptibility and potentially the thermal conductivity show overall consistency with the experimentally observed values essentially in terms of the 1D QSL. However, the signature of the 2D properties is found in a prominent peak of the spin structure factor at (π, 0) and in the structure of the ground-state wavefunction itself. In addition, the exponent of the algebraic correlation is different from the known 1D QSL found in the Heisenberg or Hubbard model. The temperature dependence of NMR T À1 1 at T < 1 K (∝T 2 ) supports vanishing density of states of spin excitation around zero energy and implies the existence of highly anisotropic but point-like nodes consistently with the present results. The present QSL offers a unified picture that bridges the 1D and 2D QSL.
It is desired to calculate the dynamical spin structure factor in the future to more directly clarify the spin dynamics with nodal excitation suggested by the structure of the wavefunction. An intriguing issue is the relation to recent studies on the charge dynamics above 2K for β 0 -EtMe 3 Sb[Pd(dmit) 2 ] 2 37 . It might be associated with essentially the 1D-like spin liquid, where the dynamical singlet formation may couple to the dimerization fluctuation of the lattice. The spin transport can be measured by attaching ferromagnetic metal leads and by estimating the difference of the spin conduction between the cases of parallel and anti-parallel magnetization for the anode and cathode. The ac response may also help to estimate the spin Drude weight studied in the present study without the sensitivity to the misalignment or randomness. These are left for future studies to stringently verify the relevance of the present finding in the collaboration of experiment and theory.

METHODS Effective Hamiltonian parameters
We have used the parameters of the effective Hamiltonian (1) derived for the low-temperature experimental crystal structure below 10 K and construct Hamiltonians on a 2D plane. Here the parameters are listed in Table 1 for the self-contained description. The interaction is screened by the interlayer screening, which must be taken into account when one solves a single layer Hamiltonian. It is established that this dimensional downfolding effect is well represented by a constant reduction of all the interactions by the amount Δ DDF from the values for the three dimensional Hamiltonian, if the subtracted value exceeds zero, and by truncation to zero if the subtracted value becomes negative 16 . For dmit compounds, it was estimated as Δ DDF = 0.18 eV. The interaction beyond the 3rd neighbor site (namely beyond V 4 and its geometrically equivalent sites) becomes small after the dimensional downfolding and we ignore them. We take into account up to the first neighbor transfer in eachã;b andc direction. The exchange interaction is also ignored because it is not expected to play any visible role in the present study. In this work, we analyze the above Hamiltonian in the form of the lattice structure shown in Fig. 6 at half filling in the canonical ensemble, with N s = L × L sites. We consider systems under the antiperiodicperiodic boundary condition. In all the cases of the compounds studied, the ground state is the Mott insulator.

Variational Monte Carlo method
In our simulations, we used a variational Monte Carlo method [21][22][23] . Our variational wavefunction takes the following form: Here, Jc ij n i n j Þ and P Js ¼ expð P i;j α Js ij S z i S z j Þ are the Gutzwiller factor 38 , the long-range Jastrow correlation factors 39,40 , and the long-range spin Jastrow correlation factor 41 , respectively. Here, S z i ¼ n i" À n i# . In practice, we impose the translational symmetry on the Gutzwiller and Jastrow factors.
To enhance the accuracy and to make the variance extrapolation explained in the next subsection easy, we also used two types of restricted Boltzmann machine (RBM) correlators 23,42,43 M s and M c , which are defined in the following equations: Here, t in Eq. (5) represents a type of the RBM correlators, i.e. t = c or s. N t α denotes the ratio of the number of the hidden units in the hidden layer to the number of the physical sites N s . x j i is a real space configuration of electrons in the sector where the total S z is zero. a, b, and W are complex variational parameters. For the measurements of the physical quantities defined in Methods, we use N c α ¼ 4 and N c α ¼ N s α ¼ 2 for the nonmagnetic states and the antiferromagnetic state, respectively. ϕ j i in Eq. (3) is a fermionic wavefunction. As ϕ j i, the generalized pairing wavefunction is employed, which is defined by where f iσ;jσ 0 are variational parameters and N e is the total number of electrons. The spin indices σ and σ 0 can be arbitrary with singlet and triplet combinations. This can accommodate the Hartree-Fock-Bogoliubov type wavefunction with antiferromagnetic (AF), charge and superconducting orders 21,44 , and flexibly describes these states and further paramagnetic metals as well. In this study, we extend it by introducing the dependence on the local density of x j i. Our extended pairing wavefunction is defined t n and V n represent the hopping parameters and Coulomb interactions, respectively, defined in Fig. 6. The units of t n , U, and V n are meV. After the dimensional downfolding, all the interaction parameters U and V n should be reduced with the amount of 180 meV.
K. Ido et al.
as follows: M nm ðxÞ ¼ f ext rnσ;rmσ 0 ðxÞ À f ext rmσ 0 ;rnσ ðxÞ; (10) where n and m are electron's indices in the sample x j i. PfM is the Pfaffian of a skew-symmetric matrix M. This extension improves accuracy of the fermionic part of the trial wavefunction especially for nonmagnetic states. We treated f ss iσ;jσ 0 , f sd iσ;jσ 0 , f ds iσ;jσ 0 , and f dd iσ;jσ 0 as complex variational parameters. In order to reduce the computational cost by saving the number of independent variational parameters, we assume that f ext ij have a sublattice structure such that f ext ij depends on the relative vector r i − r j and a sublattice index of the site j which is denoted as η(j). Thus, we can rewrite it as f ext ηðjÞ ðr i À r j Þ. In the present study on the nonmagnetic gapless states, we assumed a fully translational invariance (1 × 1 sublattice structure) and do not optimize triplet pairings (represented by the caseσ ¼ σ 0 ). For studies on the AF states, we extended the sublattice structure of f ext ij to 2 × 2. We did not use P Js and M s for the optimization of the nonmagnetic states. All the variational parameters are simultaneously optimized by using the stochastic reconfiguration method 45 .

Variance extrapolation
The AF and quantum spin liquid (QSL) states as well as the valence bond state are severely competing. To determine the lowest energy state among them, we performed extrapolations of energies to the limit of the zero energy variance [45][46][47] . For this purpose, we combined with the restricted Boltzmann machine algorithm 23,42 together with the 1st Lanczos step 48 . To obtain the ground-state energy in the zero-variance limit, we perform the linear regression by using N c α ¼ 2 and 4 for nonmagnetic states. For the collinear antiferromagnetic state, we use N c α ¼ N s α ¼ 0 and 2. Examples of the extrapolations in the present studies are found in Fig. 7, where the variance dependence of the ground-state energy for the ab initio Hamiltonian for X=EtMe 3 Sb is plotted. The limit of the zero variance is our estimate of the ground-state energy. We find severe competition between the AF and the QSL. The energies of the correlated metal and the dimer state with a finite spin gap are much higher than those of the AF and the QSL. We note that the error bars in Fig. 2 of the main text are obtained based on the propagation of errors as ε ¼ , where ε AF and ε QSL are the error arising from the linear regression for the AF and QSL states, respectively.

Physical quantities
To elucidate the nature of the ground state, we analyze the following quantities: The primarily important quantity is the spin correlation  Table 1. b Definitions of deformed triangular lattice and x − y axes. c: Definitions of screened Coulomb interactions on the lattice corresponding to non-deformed lattice in (a). U and V n (n = 1-9) on circles denote the screened Coulomb interactions from the reference site (square in the center). In the actual calculation, the same deformed structure as b is used. The values are listed in Table 1. We take into account all the transfers and interactions illustrated here.
between the site r and r þ r 0 and its Fourier transform called the spin structure factor where S r is the spin operator at the site r.
The spin Drude weight in u direction (u = x, y) is defined by the second derivative of the energy with respect to the vector potential 49 where A is inserted to the transfer by replacing as t ijσ ! t ijσ exp½iσπA Á ðr i À r j Þ=L: The spin Drude weight is associated with the spin conductivity σ(ω) as with the cutoff Λ to represent the weight around ω = 0. If the peak around ω = 0 is given by the delta function, it is reduced to We compute the statistical error of a physical quantity O arising from the Monte Carlo sampling estimated as where O h i i is the expectation value of O at the ith bin in the Monte Carlo sampling. In general, we use about 10 6 -10 7 Monte Carlo samples in each bin. N bin is the total number of bins and we typically set N bin = 5-10.

Gapless structure of wavefunction
The method to analyze the structure of Fourier transform of the singlet pairing amplitude, f ij = f i↑,j↓ , denoted by f(k) is discussed here. Although the correlation factors C in Eq. (4) and the density-dependent pairing amplitude in Eq. (9) largely modify the wavefunction character, the nodal structure is expected to be governed by f ij in Eq. (8) To obtain the fitted f(k) defined in Eq. (20), we use the following form of ϵðkÞ andΔðkÞ: ΔðkÞ ¼ 2Δ a cos k x þΔ b cosðk x þ k y Þ þΔ c cos k y À Á : We also introduce the uniform scale factor C of f(k) because the original pairing amplitudes have the ambiguity for its norm. The parameters t b ;t c ;μ;Δ a ;Δ b ,Δ c and C are simultaneously optimized by using the differential evolution method 51 implemented in SciPy 52 , which is one of the global optimization methods. We confirmed that the optimized results are similar even when we use other optimization methods such as the simplicial homology global optimization 53 . We ignore the pairing amplitudes smaller than 10 −5 during the optimization. The fitting bỹ ϵðkÞ andΔðkÞ reproduces the original data Fig. 8b quite well and is not distinguishable in the color plot. Note that the original data are obtained from the optimized P G P Jc ϕ pair for σ ≠ σ 0 with the real variational parameters α G i ; α Jc ij and f ij . Figure 8b shows the fittedεðkÞ andΔðkÞ by using Eq. (20). Since the excitation of the Hamiltonian is represented by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffif ϵðkÞ 2 þΔðkÞ 2 q , the gapless point appears at the cross points ofεðkÞ ¼ 0 andΔðkÞ ¼ 0, as are shown as red circles. Though change in quantitative slopes of dispersion and broadening may take place, these nodal points may not be seriously  Fig. 7 Example of the variance extrapolation of the energy to determine the ground state. The example is shown for the ab initio Hamiltonian for X=EtMe 3 Sb. The variance is defined as e represents the energy per site. Red circles, blue squares, green triangles and purple diamonds represent the energies of the QSL, the AF, a metal, and a spin-gap state with a finite dimer order along x direction, respectively. Open and filled symbols are obtained by using the VMC method and the VMC combined with 1st-step Lanczos method, respectively, to each of which we also apply the restricted Boltzmann machine, which enables the lower energy. Dashed lines are fitted lines for the variance extrapolations. We impose the 1 × 1 and 2 × 1 sublattice structures on the trial wavefunction for the metal and the dimer state, respectively. In all the data plot, error bars obtained in the VMC calculation are smaller than the symbol size. altered by the correlation factors such as the Gutzwiller factor P G30 . The Fourier transform of f ij denoted by f(k) can be associated with the excitation spectra 50 through the fitting to the ground state of a mean-field BCS Hamiltonian with the d-wave type superconducting order. Since the quasiparticle excitation of the BCS Hamiltonian corresponds to the spin-1/ 2 spinon, the excitation of the QSL is inferred to be characterized by the spin-1/2 Dirac-type spinon around the nodes of the d-wave superconducting state at around (±π/2, 0) and (±π/2, π). The spinon is an excitation resulted from the fractionalization of the spin and is confined. Measurable ordinary triplet excitations are given by a combination of two spinons generating the gapless points for the triplet excitations at around (0, 0), (±π, 0), (0, ±π), and (±π, ±π).

Strong coupling picture
The effective exchange coupling in the strong coupling expansion in terms of either t a /(U − V 1 ), t b /(U − V 2 ) or t c /(U − V 3 ) can be easily derived using J a ¼ 4t 2 a =ðU À V 1 Þ etc., which is summarized in Table 2. For X=EtMe 3 Sb, J a = 30.5 meV, J b = 19.6 meV, and J c = 14.4 meV (we ignore the contribution from the direct exchange because its small and more or less isotropic values less than 3 meV 18 , which have only minor effect): The onedimensionality is in fact as a consequence of not J b /J a (or J c /J a ) but small (J c − J b )/J a . Nonetheless, this factor is 0.17, which is substantially larger than the value by Kenny et al., 0.053 36 . In fact, with our present ab initio parameters, their way of estimate of the phase diagram using the above J a , J b , J c and their estimate of the ring exchange interaction K 4 ¼ 80t 2 b t 2 c =ðU À ðV 1 þ V 2 þ V 3 Þ=3Þ 3 giving K 4 /J a~0 .11 indicate that the EtMe 3 Sb[Pd(dmit) 2 ] 2 is located near the border of the antiferromagnetic phase but in the QSL side, consistently with our full quantum calculation. The estimates of the ground states for the other compounds based on their criterion are also consistent with our results. However, their recent estimate of the instability to the antiferromagnetic order contradicts our results 24 . It could be related to the neglect of the ring exchange interaction and partly because of the limitation of the random phase approximation based on the quasi-one-dimensionality employed by them 24 . In any case, it is clear that the present analysis is quantitatively the most comprehensive and accurate analysis because all orders of expansion even beyond the ring exchange from the viewpoint of the strong coupling expansion are included and the present original itinerant-electron treatment is the most strict first principles analysis.

DATA AVAILABILITY
Correspondence and requests for materials should be addressed to imada@g.ecc.utokyo.ac.jp or ido@issp.u-tokyo.ac.jp.

CODE AVAILABILITY
Code used in the present paper has been developed based on an open source software at https://www.pasums.issp.u-tokyo.ac.jp/mvmc/en/ and the algorithm is described in ref. 22 . Corrrespondence and requests for codes should be addressed to imada@g.ecc.u-tokyo.ac.jp or ido@issp.u-tokyo.ac.jp.  J n representes the effective coupling between the nearest neighbor sites along n direction. The unit of J n is meV. J b /J a , J c /J a , and (J c − J b )/J a are also listed.