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

Organic solids offer plenty of playgrounds of conspicuous phenomena including superconductivity competing with charge and spin orders, and semi-metallic excitations with Dirac dispersions RevModPhys.89.025003 ; kato2014_BCSJ . 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, $\kappa$-(ET)${}_{2}X$ (ET= bis (ethylenedithio) tetrathiafulvalene) with the anion $X$=Cu₂(CN)₃PhysRevLett.91.107001 and $\beta^{\prime}$-$X$[Pd(dmit)₂]₂ (dmit=1,3-dithiole-2-thione-4,5-dithiolate) with the cation $X$=EtMe₃Sb Itou2008 , 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)₂ 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 RevModPhys.89.025003 ; kanoda2011 ; kato2014_BCSJ . 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)₂]₂. In fact, it shows Itou2008 ; Itou2010 smaller broadening of NMR spctra than $\kappa$-(ET)₂Cu₂(CN)₃RevModPhys.89.025003 ; PhysRevB.73.140407 implying smaller effects of extrinsic spatial inhomogeneity. At low temperatures of the QSL material, EtMe₃Sb[Pd(dmit)₂]₂, the specific heat Yamashita2011 and the thermal conductivity $\kappa$ Yamashita2010 are reported to be proportional to temperature $T$, although controversies exist for the latter Yamashita2020 ; PhysRevX.9.041051 ; PhysRevLett.123.247204 . The magnetic susceptibility stays at a nonzero constant Watanabe2012 . 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 the nuclear magnetic resonance (NMR), $1/T_{1}$ seems to be scaled by $T^{2}$ in the range 0.05K$\leq T<1$K Itou2010 in contrast to the conventional Fermi liquid behavior $\propto 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 $\beta^{\prime}$-$X$[Pd(dmit)₂]₂ with $X$=Me₄P, Me₄As, Me₄Sb, Et₂Me₂As exhibiting the AF order and with $X$=EtMe₃Sb exhibiting QSL phases. We do not consider $\beta^{\prime}$-Et₂Me₂Sb[Pd(dmit)₂]₂, because the observed charge order with the simultaneous spontaneous lattice distortion in this compound nakao2005structural ; tamura2005spectroscopic 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.

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 Nakamura2012 (see Methods). The maximally localized Wannier function of the HOMO orbital, whose band crosses the Fermi level is constructed for the molecular orbitalmarzari1997 ; souza2001 . 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)₂ dimers, where a Wannier orbital is extended over a dimer.

The resultant ab initio single-band effective Hamiltonian derived in Ref. PhysRevResearch.2.032072 ; Misawa2021 has the form

	$\displaystyle{\cal H}=-\sum_{\sigma}\sum_{i,j}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}+\sum_{i}Un_{i\uparrow}n_{i\downarrow}$
	$\displaystyle+\sum_{i<j}V_{ij}(n_{i}-1)(n_{j}-1),$		(1)

where $i$, $j$ represent the dimer indices, and $c_{i\sigma}^{\dagger}$ ($c_{i\sigma}$) is the creation (annihilation) operator of electrons with spin $\sigma$ (=$\uparrow$ or $\downarrow$) at the $i$-th Wannier orbital, and the number operator is $n_{i}=\sum_{\sigma}n_{i\sigma}$ with $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$. Here, $t_{ij}$ is the hopping parameters depending on the relative coordinate vector $\bm{r}_{i}-\bm{r}_{j}$, where $\bm{r}_{i}$ is the position vector of the center of the $i$-th Wannier orbital. In the present study, for $t_{ij}$ and $V_{ij}$, we retain the nearest neighbor pair of $\bm{r}_{i}$ and $\bm{r}_{j}$ in each direction of $\tilde{a}$, $\tilde{b}$ and $\tilde{c}$ as is illustrated in Fig. 1. 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 Tahara2008 ; MISAWA2019447 ; Nomura2017 to the ab initio Hamiltonians to reach highly accurate ground states. For details of the numerical method, see Methods.

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 $\Delta e$ shows that $(t_{c}-t_{b})/t_{a}$ is indeed a principally important parameter for 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 projected BCS study of the Hubbard-type model by using the ab initio hopping PhysRevB.88.155139 .

The spin-spin correlations in real space $C(\bm{r})$ for the QSL state of $X$=EtMe₃Sb shown in Figure 3(d) with the log-log plot in the $x$ (namely, $t_{a}$) direction indicates a power law $C(\bm{r})\propto|\bm{r}|^{-p}$ with $p=1.88\pm 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({\bm{r}})\propto\exp[-|{\bm{r}}|/\xi_{\perp}]$ with the correlation length $\xi_{\perp}\sim 1.33\pm 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}$ (namely $\tilde{b}$) and $t_{c}$ (namely, $\tilde{c}$) directions. The anisotropic correlation develops the ridge in the $q_{y}$ dependence of $S(\bm{q})$.

Finally, the spin Drude weight $D_{\rm s}$ calculated in the QSL state for $X$=EtMe₃Sb is shown in Fig. 4. The size dependence of $D_{\rm s}$ for the QSL shown in Fig. 4(b) assures that $D_{\rm 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₃Sb, which is our second central result of the paper.

One might regard the present QSL as essentially the same as the 1D spin liquid PhysRevB.89.235107 ; PhysRevB.74.012407 ; PhysRevB.74.014408 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 $(\pi,0)$ in $S(\bm{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 $(\pi,0)$-$(\pi,\pi)$ 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 $(\pi,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(\bm{q})$ toward $(\pi,\pi)$ from $(\pi,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 $(\pm\pi/2,\pm\pi/2)$ Ferrari_2018 ; NomuraJ1J2_2021 . Let us discuss the possibility of the fractionalization for the present QSL of $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂. 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 $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂ 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 $(\pm\pi/2,0)$ and $(\pm\pi/2,\pi)$ 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)$, $(\pm\pi,0)$, $(0,\pm\pi)$, and $(\pm\pi,\pm\pi)$. This structure has an essential similarity to the 2D $J_{1}$-$J_{2}$ Heisenberg model NomuraJ1J2_2021 , 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 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

$\displaystyle H=J\sum_{i}[{S}_{i}^{x}{S}_{i+1}^{x}+{S}_{i}^{y}{S}_{i+1}^{y}+{S}_{i}^{z}{S}_{i+1}^{z}]$

(2)

with $J>0$ for the spin half operator $\bm{S}_{i}=({S}_{i}^{x},{S}_{i}^{y},{S}_{i}^{z})$ at site $i$. The specific heat of the 1D Heisenberg model is given by the $T$-linear coefficient $\gamma=2k_{\rm B}^{2}N_{\rm A}/3J$ above 1K Takahashi1973 , where $k_{\rm B}$ is Boltzmann constant and $N_{\rm 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 $\tilde{a}$ (namely, $t_{a}$) direction is estimated as $J_{a}=4t_{a}^{2}/(U-V_{1})\sim 0.031$ eV for $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂, by using the ab initio parameters listed in Table 1 of Methods. Then $\gamma$ is estimated as $\gamma\sim 15$ mJ K^-2 mole^-1, comparable with the experimental value for $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂, $\gamma^{\rm exp}\sim 15-20$ mJ K^-2mole^-1 Yamashita2011 ; PhysRevLett.123.247204 .

The uniform magnetic susceptibility of the 1D Heisenberg model is $\chi_{\rm 1D}=g^{2}\mu_{\rm B}^{2}/(J_{a}\pi^{2})$ with the $g$-factor and the Bohr magneton $\mu_{\rm B}$ Griffiths1964 ; Yang1966 . By using $J_{a}\sim 0.031$eV, and $g\sim 2$, we obtain $\chi/\mu_{\rm B}^{2}\sim 13$ eV^-1, while the experimental value is $\chi/\mu_{\rm B}^{2}\sim 12-24$ eV^-1 per the dimer Itou2008 ; Watanabe2012 , which is consistent with each other by considering the experimental uncertainty.

Frequency dependent thermal conductivity $\kappa$ of the 1D Heisenberg model is estimated to be $\kappa^{\rm 1D}(\omega)=\kappa_{\rm s}\delta(\omega)$ with $\kappa_{\rm s}/{T}={\pi^{3}J}/(6a\hbar^{2})$, where $a$ is the lattice constantKlumper2002 . If the delta function is replaced by the Lorentzian $\omega_{\rm W}/\pi(\omega^{2}+\omega_{\rm W}^{2})$ to take into account the lifetime $\tau=2\pi/\omega_{\rm W}$, we obtain $\kappa_{\rm s}^{\rm 1D}(\omega=0)/T=\pi^{2}J_{a}/(6a\hbar^{2}\omega_{\rm W})$. By considering $a\sim 10^{-9}$m and $J_{a}\sim 0.031$eV, $\kappa_{\rm s}^{\rm 1D}(\omega=0)/T$ is estimated as $7.7\times 10^{22}/\hbar\omega_{\rm W}$[Jms]^-1. Then if the experimental value reported in Ref. Yamashita2010 is employed, and the one-to-one correspondence between the 1D Heisenberg and the experimental value is assumed, $\kappa_{\rm exp}^{\rm dmit}/T\sim 0.2/k_{\rm B}^{2}$[Jms]^-1 corresponds to the carrier relaxation time $\tau=2\pi/\omega_{W}\sim 9.0\times 10^{-12}$sec. A simple expectation of the spin velocity $v_{\rm s}=J_{a}a/\hbar$ results in the mean free path $\ell_{\rm s}=v_{\rm s}\tau=0.42\mu$m, which is a value comparable to the reported estimate $\sim 1\mu$m Yamashita2010 . More precisely, the experimental value interpreted by the 1D Heisenberg model suggests $\hbar\omega_{\rm W}/J\sim 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 $\ell_{\rm s}$ is substantially longer than that of the inorganic compound Cs₂CuCl₄PhysRevResearch.1.032022 , where a 1D QSL was claimed. The present result implies that the 1D QSL found here has potentially much longer $\tau$. 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}\propto T^{2}$ below 1K Itou2010 , 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 $\sim(\pm\pi/2,0)$ and $\sim(\pm\pi/2,\pi)$, 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. $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂ is located in this parameter space close to the 1D QSL. Furthermore, the controversy and sensitivity of the thermal conductivity Yamashita2020 ; PhysRevX.9.041051 ; PhysRevLett.123.247204 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 $\kappa$ and serious sample or measurement dependence. Effects of randomness may also seriously disturb the ideal behavior of the 1D-like spin liquid PhysRevX.9.041051 . 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. PhysRevMaterials.4.044403 . (See also Methods Unconventional dual 1D-2D quantum spin liquid revealed by ab initio studies on organic solids family)

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₄P, Me₄As, Me₄Sb, Et₂Me₂As and the QSL for EtMe₃Sb. The relative stabilities of the four AF 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 $(\pi,0)$ and in the structure of the ground-state wave function 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 ($\propto 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 wave function. An intriguing issue is the relation to recent studies on the charge dynamics above 2K for $\beta^{\prime}$-EtMe₃Sb[Pd(dmit)₂]₂ PhysRevB.97.035131 . 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.

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 $\Delta_{\rm 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 Nakamura2012 . For dmit compounds, it was estimated as $\Delta_{\rm DDF}=0.18{\rm eV}$. The interaction beyond the neighboring site (namely $V_{4}$ and beyond) becomes small after the dimensional downfolding and we ignore them. We then take into account up to the first neighbor transfer and interaction in each $\tilde{a},\tilde{b}$ and $\tilde{c}$ 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=L\times L$ sites. We consider systems under the antiperiodic-periodic boundary condition. In all the cases of the compounds studied, the ground state is the Mott insulator.

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) correlatorsCarleo2017 ; Nomura2017 ; Ferrari2019 ${\cal M}^{\rm s}$ and ${\cal M}^{\rm c}$, which are defined in the following equations:

$${\cal M}^{t}=e^{a\sum_{i,\sigma}p^{t}_{i,\sigma}}\prod_{\alpha}^{N^{t}_{\alpha}}\prod_{r}^{N}2\cosh\left(b_{\alpha}+\sum_{i,\sigma}W_{i+r,\sigma,\alpha}p^{t}_{i,\sigma}\right),$$

(5)

$$p^{s}_{i,\sigma}=S_{i}^{z}=\begin{cases}1&\text{for $S_{i}^{z}\mathinner{|{x}\rangle}=\mathinner{|{x}\rangle}$}\\ 0&\text{for $S_{i}^{z}\mathinner{|{x}\rangle}=0$}\\ -1&\text{for $S_{i}^{z}\mathinner{|{x}\rangle}=-\mathinner{|{x}\rangle}$}\end{cases},$$

(6)

$$p^{c}_{i,\sigma}=2(n_{i}-1)^{2}-1=\begin{cases}1&\text{for $n_{i}\mathinner{|{x}\rangle}\neq\mathinner{|{x}\rangle}$}\\ -1&\text{for $S_{i}^{z}\mathinner{|{x}\rangle}\neq 0$}\end{cases}.$$

(7)

Here, $t$ in Eq. (5) represents a type of the RBM correlators, i.e. $t=c$ or $s$. $N^{t}_{\alpha}$ denotes the ratio of the number of the hidden units in the hidden layer to the number of the physical sites $N$. $\mathinner{|{x}\rangle}$ is a real space configuration of electrons in the sector where the total $S^{z}$ is zero. $a$, $\bm{b}$, and $\bm{W}$ are complex variational parameters. For the measurements of the physical quantities defined in Methods Unconventional dual 1D-2D quantum spin liquid revealed by ab initio studies on organic solids family, we use $N^{c}_{\alpha}=4$ and $N^{c}_{\alpha}=N^{s}_{\alpha}=2$ for the nonmagnetic states and the antiferromagnetic state, respectively.

$|\phi\rangle$ in Eq. (3) is a fermionic wavefunction. As $|\phi\rangle$, the generalized pairing wave function is employed, which is defined by

$\displaystyle|\phi^{\rm pair}\rangle=\left(\sum_{i\sigma,j\sigma^{\prime}}f_{i\sigma,j\sigma^{\prime}}c_{i\sigma}^{\dagger}c_{j\sigma^{\prime}}^{\dagger}\right)^{N_{\rm e}/2}|0\rangle,$

(8)

where $f_{i\sigma,j\sigma^{\prime}}$ are variational parameters and $N_{\rm e}$ is the total number of electrons. The spin indices $\sigma$ and $\sigma^{\prime}$ can be arbitrary with singlet and triplet combinations. This can accommodate the Hartree-Fock-Bogoliubov type wave function with antiferromagnetic (AF), charge and superconducting ordersTahara2008 ; giamarchi1991 , 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 $\mathinner{|{x}\rangle}$. Our extended pairing wavefunction is defined as follows:

	$\displaystyle\mathinner{|{\phi^{\rm ext}}\rangle}=\sum_{x}{\rm Pf}M(x)\mathinner{|{x}\rangle},$		(9)
	$\displaystyle M_{nm}(x)=f^{\rm ext}_{\bm{r}_{n}\sigma,\bm{r}_{m}\sigma^{\prime}}(x)-f^{\rm ext}_{\bm{r}_{m}\sigma^{\prime},\bm{r}_{n}\sigma}(x),$		(10)

$$f^{\rm ext}_{i\sigma,j\sigma^{\prime}}(x)=\begin{cases}f_{i\sigma,j\sigma^{\prime}}^{\rm ss}&\text{for $n_{i,-\sigma}\mathinner{|{x}\rangle}=0$, $n_{j,-\sigma^{\prime}}\mathinner{|{x}\rangle}=0$}\\ f_{i\sigma,j\sigma^{\prime}}^{\rm sd}&\text{for $n_{i,-\sigma}\mathinner{|{x}\rangle}=0$, $n_{j,-\sigma^{\prime}}\mathinner{|{x}\rangle}=\mathinner{|{x}\rangle}$}\\ f_{i\sigma,j\sigma^{\prime}}^{\rm ds}&\text{for $n_{i,-\sigma}\mathinner{|{x}\rangle}=\mathinner{|{x}\rangle}$, $n_{j,-\sigma^{\prime}}\mathinner{|{x}\rangle}=0$}\\ f_{i\sigma,j\sigma^{\prime}}^{\rm dd}&\text{for $n_{i,-\sigma}\mathinner{|{x}\rangle}=\mathinner{|{x}\rangle}$, $n_{j,-\sigma^{\prime}}\mathinner{|{x}\rangle}=\mathinner{|{x}\rangle}$}\end{cases},$$

(11)

where $n$ and $m$ are electron’s indices in the sample $\mathinner{|{x}\rangle}$. ${\rm Pf}M$ 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_{i\sigma,j\sigma^{\prime}}^{\rm ss}$, $f_{i\sigma,j\sigma^{\prime}}^{\rm sd}$, $f_{i\sigma,j\sigma^{\prime}}^{\rm ds}$, and $f_{i\sigma,j\sigma^{\prime}}^{\rm dd}$ as complex variational parameters.

In order to reduce the computational cost by saving the number of independent variational parameters, we assume that $f^{\rm ext}_{ij}$ have a sublattice structure such that $f^{\rm ext}_{ij}$ depends on the relative vector ${\bm{r}}_{i}-{\bm{r}}_{j}$ and a sublattice index of the site $j$ which is denoted as $\eta(j)$. Thus, we can rewrite it as $f^{\rm ext}_{\eta(j)}({\bm{r}}_{i}-{\bm{r}}_{j})$. In the present study on the nonmagnetic gapless states, we assumed a fully translational invariance (1$\times$1 sublattice structure) and do not optimize triplet pairings (represented by the case$\sigma=\sigma^{\prime}$). For studies on the AF states, we extended the sublattice structure of $f^{\rm ext}_{ij}$ to $2\times 2$. We did not use ${\cal P}^{\rm J_{s}}$ and ${\cal M}^{\rm s}$ for the optimization of the nonmagnetic states. All the variational parameters are simultaneously optimized by using the stochastic reconfiguration methodsorella2001 .

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 kwon1993 ; imada2000 ; sorella2001 . For this purpose, we combined with the restricted Boltzmann machine algorithm Carleo2017 ; Nomura2017 together with the 1st Lanczos stepheeb1993systematic . To obtain the ground-state energy in the zero-variance limit, we perform the linear regression by using $N^{c}_{\alpha}=2$ and $4$ for nonmagnetic states. For the collinear antiferromagnetic state, we use $N^{c}_{\alpha}=N^{s}_{\alpha}=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₃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 $\varepsilon=\sqrt{\varepsilon_{\rm AF}^{2}+\varepsilon_{\rm QSL}^{2}}$, where $\varepsilon_{\rm AF}$ and $\varepsilon_{\rm QSL}$ are the error arising from the linear regression for the AF and QSL states, respectively.

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 Kohn1964

$\displaystyle D^{u}_{\rm s}=\frac{1}{4\pi}\frac{\partial^{2}E}{\partial A_{u}^{2}},$

(14)

where $\bm{A}$ is inserted to the transfer by replacing as

$\displaystyle t_{ij\sigma}\rightarrow t_{ij\sigma}\exp[i\sigma\pi\bm{A}\cdot(\bm{r}_{i}-\bm{r}_{j})/L].$

(15)

The spin Drude weight is associated with the spin conductivity $\sigma(\omega)$ as

$\displaystyle D_{\rm s}=\int_{0}^{\Lambda}\sigma(\omega)$

(16)

with the cutoff $\Lambda$ to represent the weight around $\omega=0$. If the peak around $\omega=0$ is given by the delta function, it is reduced to

$\displaystyle\sigma(\omega)=D_{\rm s}\delta(\omega).$

(17)

We compute the statistical error of a physical quantity $O$ arising from the Monte Carlo sampling estimated as

$\displaystyle\sigma_{\rm err}(O)=\frac{1}{\sqrt{N_{\rm bin}}}\sqrt{\frac{\sum_{i}[\mathinner{\langle{O}\rangle}_{i}-\frac{1}{N_{\rm bin}}\left(\sum_{i}\mathinner{\langle{O}\rangle}_{i}\right)]^{2}}{N_{\rm bin}-1}},$

where $\mathinner{\langle{O}\rangle}_{i}$ is the expectation value of $O$ at the $i$-th bin in the Monte Carlo sampling. In general, we use about $10^{6}$-$10^{7}$ Monte Carlo samples in each bin. $N_{\rm bin}$ is the total number of bins and we typically set $N_{\rm bin}=5$-$10$.

The method to analyze the structure of Fourier transform of the singlet pairing amplitude, $f_{ij}=f_{i\uparrow,j\downarrow}$, denoted by $f(\bm{k})$ is discussed here. Although the correlation factors ${\cal 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) for $\sigma\neq\sigma^{\prime}$. The Fourier transform $f(\bm{k})$ can be interpreted as the solution of the BCS mean-field Hamiltonian

	$\displaystyle{\cal H}_{\rm BCS}=-\sum_{\bm{k},\sigma}\tilde{\epsilon}(\bm{k})c_{\bm{k}\sigma}^{\dagger}c_{\bm{k}\sigma}+\sum_{k}\tilde{\Delta}(\bm{k})[c_{\bm{k}\sigma}^{\dagger}c_{-\bm{k}-\sigma}^{\dagger}+H.c]$
			(19)

as PhysRevLett.87.097201

$\displaystyle f(\bm{k})=\frac{\tilde{\Delta}(\bm{k})}{\tilde{\epsilon}(\bm{k})+\sqrt{\tilde{\epsilon}(\bm{k})^{2}+\tilde{\Delta}(\bm{k})^{2}}}.$

(20)

To obtain the fitted $f(\bm{k})$ defined in Eq. (20), we use the following form of $\tilde{\epsilon}(\bm{k})$ and $\tilde{\Delta}(\bm{k})$:

	$\displaystyle\tilde{\epsilon}(\bm{k})=-2\left(\cos k_{x}+\tilde{t_{b}}\cos(k_{x}+k_{y})+\tilde{t_{c}}\cos k_{y}\right)-\tilde{\mu},$		(21)
	$\displaystyle\tilde{\Delta}(\bm{k})=2\left(\tilde{\Delta_{a}}\cos k_{x}+\tilde{\Delta_{b}}\cos(k_{x}+k_{y})+\tilde{\Delta_{c}}\cos k_{y}\right).$		(22)

We also introduce the uniform scale factor $C$ of $f(\bm{k})$ because the original pairing amplitudes have the ambiguity for its norm. The parameters $\tilde{t_{b}},\tilde{t_{c}},\tilde{\mu},\tilde{\Delta_{a}},\tilde{\Delta_{b}}$, $\tilde{\Delta_{c}}$ and $C$ are simultaneously optimized by using the differential evolution methodStorn1997 implemented in SciPy2020SciPy-NMeth , 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 optimizationEndres2018 . We ignore the pairing amplitudes smaller than $10^{-5}$ during the optimization. The fitting by $\tilde{\epsilon}(\bm{k})$ and $\tilde{\Delta}(\bm{k})$ reproduces the original data Fig. 8(b) quite well and is not distinguishable in the color plot. Note that the original data are obtained from the optimized ${\cal P^{\rm G}}{\cal P^{\rm J_{c}}}\mathinner{|{\phi^{\rm pair}}\rangle}$ for $\sigma\neq\sigma^{\prime}$ with the real variational parameters $\alpha^{\rm G}_{i},\alpha^{\rm J_{c}}_{ij}$ and $f_{ij}$.

Figure 8(b) shows the fitted $\tilde{\epsilon}(\bm{k})$ and $\tilde{\Delta}(\bm{k})$ by using Eq. (20). Since the excitation of the Hamiltonian is represented by $\sqrt{\tilde{\epsilon}(\bm{k})^{2}+\tilde{\Delta}(\bm{k})^{2}}$, the gapless point appears at the cross points of $\tilde{\epsilon}(\bm{k})=0$ and $\tilde{\Delta}(\bm{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 altered by the correlation factors such as the Gutzwiller factor $\cal{P}^{\rm G}$Ferrari_2018 . The Fourier transform of $f_{ij}$ denoted by $f({\bm{k}})$ can be associated with the excitation spectra PhysRevLett.87.097201 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 $(\pm\pi/2,0)$ and $(\pm\pi/2,\pi)$. 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)$, $(\pm\pi,0)$, $(0,\pm\pi)$, and $(\pm\pi,\pm\pi)$.

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_{a}^{2}/(U-V_{1})$ etc., which is summarized in Table 2. For $X$=EtMe₃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 3meV Misawa2021 , which have only minor effect): The one-dimensionality 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 PhysRevMaterials.4.044403 . 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_{b}^{2}t_{c}^{2}/(U-(V_{1}+V_{2}+V_{3})/3)^{3}$ giving $K_{4}/J_{a}\sim 0.11$ indicate that the EtMe₃Sb[Pd(dmit)₂]₂ 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 resultsPhysRevMaterials.5.084412 . 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 PhysRevMaterials.5.084412 . 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.

The authors thank Reizo Kato for discussions and clarifications on the experimental results and Shigeki Fujiyama for notifying his unpublished result on the coexistence of AF and QSL phases for $\beta^{\prime}$-Et₂Me₂As[Pd(dmit)₂]₂ and $\beta^{\prime}$-Me₄Sb[Pd(dmit)₂]₂. KI thanks Yusuke Nomura for useful comments and the implementation on the restricted Boltzmann machine correlator. TM and KY thank Takao Tsumuraya for discussions on the derivation of the $ab$ $initio$ Hamiltonians. This work was supported in part by KAKENHI Grant No. 16H06345 and 19K14645 from JSPS. This research was also supported by MEXT as “program for Promoting Researches on the Supercomputer Fugaku”(Basic Science for Emergence and Functionality in Quantum Matter - Innovative Strongly Correlated Electron Science by Integration of Fugaku and Frontier Experiments -, JPMXP1020200104). We thank the Supercomputer Center, the Institute for Solid State Physics, The University of Tokyo for the use of the facilities. We also thank the computational resources of supercomputer Fugaku provided by the RIKEN Center for Computational Science (Project ID: hp200132 and hp210163) and Oakbridge-CX in the Information Technology Center, The University of Tokyo.