Mechanism of superconductivity and electron-hole doping asymmetry in κ-type molecular conductors

Unconventional superconductivity in molecular conductors is observed at the border of metal-insulator transitions in correlated electrons under the influence of geometrical frustration. The symmetry as well as the mechanism of the superconductivity (SC) is highly controversial. To address this issue, we theoretically explore the electronic properties of carrier-doped molecular Mott system κ-(BEDT-TTF)2X. We find significant electron-hole doping asymmetry in the phase diagram where antiferromagnetic (AF) spin order, different patterns of charge order, and SC compete with each other. Hole-doping stabilizes AF phase and promotes SC with dxy-wave symmetry, which has similarities with high-Tc cuprates. In contrast, in the electron-doped side, geometrical frustration destabilizes the AF phase and the enhanced charge correlation induces another SC with extended-s + \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_{x^2 - y^2}$$\end{document}dx2-y2wave symmetry. Our results disclose the mechanism of each phase appearing in filling-control molecular Mott systems, and elucidate how physics of different strongly-correlated electrons are connected, namely, molecular conductors and high-Tc cuprates.

U nderstanding the intimate correlation among metalinsulator (MI) transition, magnetism, and superconductivity (SC) is one of the most challenging issues in modern condensed matter physics. The most well-studied example is the high-T c cuprates, where SC is observed when mobile carriers are doped into the parent antiferromagnetic (AF) Mott insulators [1][2][3] . A general understanding there, supported by various experiments and theories, is that strong AF spin fluctuation mediates the d-wave SC that appears through the fillingcontrol Mott MI transition generating mobile charge carriers. However, can we export this mechanism to other strongly correlated materials? To address this question, it is crucial to make a comparison among different classes of materials. In this respect, heavy fermion compounds and molecular conductors provide such opportunities [4][5][6][7][8] .
The family of quasi two-dimensional molecular conductors κ-(ET) 2 X (ET = BEDT-TTF, and X takes different monovalent anions 9 ) is in fact compared often with the cuprates 10 . They indeed have common factors: simple quasi-two-dimensional electronic structure to begin with in the non-interacting limit and the Mott MI transition and SC closely related with each other. However, there are important differences: First, in κ-(ET) 2 X, SC appears through the bandwidth-control Mott transition; the carrier density is usually unchanged but the pressure (either physically or chemically) is the controlling factor. Although the variation of the carrier density is necessary for direct comparisons, it has not been realized in κ-(ET) 2 X for a long time due to experimental difficulties. Second, while the cuprates are basically governed by the physics nearby 1/2-filling, κ-(ET) 2 X is a 3/4-filled system. The similarity enters when the so-called dimer approximation is applied in the latter 11 , resulting in the effective 1/2-filled system (dimer model). Although the dimer model has been extensively studied using various theoretical methods [12][13][14][15][16][17][18][19][20][21][22][23][24][25] , the validity of the dimer approximation itself is recently reexamined [26][27][28][29][30][31][32] . Especially, the importance of the intradimer charge degree of freedom 33,34 and intersite Coulomb interactions 27,29 , which are discarded in the dimer approximation, has attracted much attention because of recent experimental suggestions [35][36][37][38] . Third, in κ-(ET) 2 X, SC is observed not only next to the AF insulators but also to nonmagnetic (candidate of gapless spin-liquid) insulators 8,[39][40][41] . The strong influence of geometrical frustration owing to the anisotropic triangular arrangement of dimers is present in this family.
Recently, carrier doping has been realized either chemically in κ-(ET) 4 Hg 3−δ Y 8 (Y = Br or Cl) [42][43][44][45] or in κ-(ET) 2 Cu[N(CN) 2 ]Cl (κ-Cl) by using electric-double-layer transistor (EDLT) technique 46,47 , revealing intriguing phenomena such as a domeshaped SC region, anomalous metallic behaviors, and significant electron-hole doping asymmetry, which are all reminiscent of the high-T c cuprates. Therefore, κ-type ET systems can now provide a unique playground of both filling-control and bandwidthcontrol Mott transitions with SC phases nearby, for which a unified theoretical understanding is highly desired.
In this paper, we theoretically study the ground-state properties of κ-(ET) 2 X varying the carrier number from 3/4-filling, in order to elucidate the electronic phases appearing near the Mott transition in this system, especially SC, and to investigate their stabilities beyond mean field treatments. The intradimer charge degree of freedom and intersite Coulomb interactions are explicitly considered. We find that the ground-state phase diagram shows significant electron-hole asymmetry in the stability of AF phase and in terms of competing two types of SC. While in the hole-doped side d xy -wave SC is favored by the AF spin fluctuation as in the high-T c cuprates, the electron doping highlights the geometrical frustration and the charge degree of freedom, which are unique in κ-(ET) 2 X, stabilizing extended-s +d x 2 Ày 2 -wave SC.
The electron-hole doping asymmetry, including the symmetry of SC, is attributed to the degree of frustration that is controlled by carrier doping. This is a conceptually new perspective to κ-(ET) 2 X, which can also be applied to other frustrated systems in general. Our results, beyond the usual description based on the 1/2-filled dimer model, thus provide new understanding of how physics of molecular conductors and high-T c cuprates are distinct.

Results
Model derivation and framework. The electronic properties of molecular conductors are modeled by a simple model where the molecules are replaced by lattice sites 48 . They are described by the extended Hubbard model (EHM) 11,49 , a textbook model for studying correlated electrons. The Hamiltonian is given as where c y iσ (c iσ ) is a creation (annihilation) operator of electron at molecular site i with spin σ(=↑,↓), n iσ ¼ c y iσ c iσ , and n i = n i↑ + n i↓ . U and V ij are on-site and intersite Coulomb repulsions, respectively. <i, j> denotes pairs of neighboring molecules in the κ-type geometry, labeled by b 1 , b 2 , p, and q, as shown in Fig. 1a.
The tight-binding parameters t ij are set for the deuterated κ-(ET) 2 Cu[N(CN) 2 ]Br (κ-Br), which locates very close to the MI transition 50 , and are adopted from a first-principles band calculation as ðt b 1 ; t b 2 ; t p ; t q Þ ¼ ð196; 65; 105; À39Þ meV ¼ ð1:0; 0:332; 0:536; À0:199Þ t b 1 51 . We set the largest hopping integral t b 1 as the unit of energy. The unit cell is a rectangle with R x × 2R y and δ ± = (δ x , ±δ y ) are vectors connecting the centers of Br. E F denotes the Fermi energy for n hole = 1/2 (undoped case). Dotted lines correspond to the Fermi energy for n hole = 1/ 3, 1/2, and 2/3. High symmetry points of momentum k are Γ(0,0), M(π/R x , π/2R y ), X(π/R x ,0), and Y(0, π/2R y ) (See also Fig. 3b) molecules facing each other in a dimer (see Fig. 1a). Here, we set 29 . The non-interacting band structure is shown in Fig. 1b. Among the four energy bands, the upper two bands (bands 1 and 2) contribute to form the Fermi surface (FS). For the undoped case, the electron density per molecular site is n = 3/2 (3/4-filling) and it corresponds to the hole density n hole = 2 − n = 1/2. In this study, we change n hole from 1/3 to 2/3 to investigate the doping dependence of the system. The corresponding Fermi energies are indicated in Fig. 1b by dotted lines. The effect of Coulomb interactions is treated using a variational Monte Carlo (VMC) method [52][53][54] . The trial wave function considered here is a Gutzwiller-Jastrow type, |Φ> is a one-body part constructed by diagonalizing the one-body Hamiltonian, and P J c and P J s are charge and spin Jastrow factors, respectively. The explicit form of them are described in Methods. In the following, we show results for 1152 molecular sites (corresponding to L = 24, see in Methods), which is large enough to avoid finite size effects.
Ground-state phase diagram. Figure 2a shows the ground-state phase diagram. The hole density n hole = 2 − n and the on-site Coulomb interaction U=t b 1 are varied as parameters, while the ratio between U and the largest intersite Coulomb interaction V b 1 is fixed at V b 1 =U ¼ 0:50. The other intersite Coulomb interactions are set as ðV b 2 ; V p ; V q Þ ¼ ð0:56; 0:66; 0:58Þ V b 1 , assuming the 1/r-dependence. At n hole = 1/2 (undoped case) 29 , a first-order phase transition occurs, with increasing U=t b 1 , from a paramagnetic metal (PM) to a dimer-type AF (DAF) phase in which the spins between dimers order in a staggered way as shown in Fig. 2b. This transition corresponds to the Mott MI transition. As U=t b 1 increases further, there appears a polar charge-ordered (PCO) phase breaking the inversion symmetry 31,34 with AF spin order, which can avoid the energy loss of V b 1 , V b 2 , and V p , at the expense of the energy loss of V q as shown schematically in Fig. 2c. The DAF and PCO phases are insulating at n hole = 1/2. They have also been found in previous studies for the 3/4-filled Hubbard models 11,31,49 and the effective strong coupling models 33,34 , and are stabilized in the relevant parameter regions for κ-(ET) 2 X. Experimentally, the DAF phase is widely observed in κ-(ET) 2 X as AF dimer-Mott insulator and the PCO phase is proposed to be related to the dielectric anomaly observed in κ-(ET) 2 Cu 2 (CN) 3 (κ-CN) 36 and the insulating phase in κ-(ET) 2 Hg(SCN) 2 Cl 55 . Note that SC is a metastable state for U=t b 1 = 7-11.5, lying on each side of the DAF-PCO boundary.
Away from n hole = 1/2, significant doping asymmetry is observed and several different phases appear. For the holedoped side (n hole > 1/2), while the PCO phase is rapidly suppressed, the DAF phase is enhanced to a smaller U=t b 1 region toward n hole = 2/3. Note that in these phases the system becomes metallic once the doping is finite. Furthermore, a 3-fold chargeordered (3-fold CO-1) phase appears for larger U=t b 1 . The 3-fold CO-1 phase shows charge disproportionation and magnetic order as shown in Fig. 2d; hole-rich sites form a two-dimensional network with AF spin order. This phase is insulating at n hole = 2/3 (along the brown line in Fig. 2a) since the electron density fits the commensurability, and metallic for other hole densities because the excess holes can move through the ordered holes.
For the electron-doped side (n hole < 1/2), the situation is much different. The PCO and DAF (both become metallic) phases are rapidly suppressed and another CO (3-fold CO-2) phase and a SC phase appear. The pattern of the charge disproportionation is opposite to that in the 3-fold CO-1 phase (hole rich ↔ hole poor) as shown in Fig. 1e; this configuration can fully avoid the intersite Coulomb interactions. Similar to the 3-fold CO-1 phase, the 3fold CO-2 phase is insulating at n hole = 1/3 (along the blue line in Fig. 2a) and metallic for other hole densities. Note that the 3-fold CO-2 phase is stabilized also for n hole = 1/2 when V b 1 =U is larger (≥0.55) 29,31 and is smoothly connected in the parameter space. While the 3-fold CO-1 phase is accompanied by the magnetic order, 3-fold CO-2 is nonmagnetic. We have tried several magnetic ordering patterns that coexist with 3-fold CO-2. However, none of them are stabilized because the CO pattern in the hole-rich sites forms a triangular-like structure and the spin degree of freedom is fully frustrated, and furthermore the distances between hole-rich sites are much longer than the original intermolecular bonds and therefore the effective magnetic exchange couplings are quite small. For smaller U=t b 1 , the SC phase is realized by doping, located between the DAF/3fold CO-2 and PM phases. The symmetry of the SC is the extended-s + d x 2 Ày 2 -wave type, same with the one shown in our previous study at n hole = 1/2 29 . Details are discussed later.
Although SC does not appear as the ground state in the holedoped side, we find finite superconducting condensation energy in the phase diagram. Figure 2f shows the region where the condensation energy is finite, ignoring other ordered phases by setting Weiss fields to be zero in |Φ〉. It is possible that the hidden SC phase appears if the DAF phase is destabilized by, e.g., disorder effect associated with doping or phase separation. Therefore, it is worthwhile to study the most favored SC phase even if it is a metastable state. While the d xy -wave SC is dominant for most of the hole-doped side, the extended-s + d x 2 Ày 2 -wave SC is stabilized for the electron-doped side. Namely, the symmetry of SC changes with carrier doping. Note that the charge correlation is greatly enhanced toward regions indicated by gray shade in Fig. 2f. In these regions, the mobility of holes are strongly restricted due to the strong intersite Coulomb interactions, and stable VMC simulations are difficult unless additional Weiss fields that induce long-range CO are introduced in |Φ〉.
Fermi surface and spin structure factor. The electron-hole doping asymmetry is closely related to the shape of the FS and the interdimer magnetic fluctuations. Figure 3a-c show the noninteracting FS for n hole = 2/3, 1/2, and 1/3. As n hole increases from 1/2 (hole doping), the FS shifts toward the right and left edges of the first Brillouin zone (M-X line in Fig. 3b). Since the energy gap of the DAF order opens along the Brillouin zone edge 11 , the DAF order becomes more favored for hole doping. For n hole = 2/3, the FS almost touches the Brillouin zone edge, and there the DAF region extends down to U=t b 1 $ 2:5, as shown in Fig. 2a. Note that the Fermi energy is located in the vicinity of van Hove singularity at n hole = 2/3 as shown in Fig. 1b. Around this hole density, anomalous behavior such as pseudogap phenomena is naively expected 46 . In clear contrast, the FS departs from the M-X line for electron doping, consistent with the tendency of the DAF order being rapidly suppressed for n hole < 1/2.
Next, Fig. 3d-f show the interdimer spin structure factor defined as, for n hole = 2/3, 1/2, and 1/3. Here, N dim (=L 2 ) is the total number of dimers and M dim l ¼ ðn 2lÀ1" þ n 2l" Þ À ðn 2lÀ1# þ n 2l# Þ is the total spin density within l-th dimer formed by molecular sites 2l − 1 and 2l with the central position r l 56 . Since the dimer centers form the anisotropic triangular lattice 11 , the corresponding first Brillouin zone is the anisotropic hexagon. As shown in Fig. 3d, S dim (q) for n hole = 2/3 peaks around (0, ±π/R y ), which corresponds to the DAF spin configuration, suggesting that the AF spin fluctuation is enhanced by the Coulomb interactions and thus the DAF order is favored. On the other hand, for n hole = 1/2, the peaks appear around six vertices of the Brillouin zone, as shown in Fig. 3e. This implies that the spin structure becomes more triangular-lattice like (frustrated) and the AF spin fluctuation is suppressed as compared with that at n hole = 2/3. For n hole = 1/3, the peak structures almost diminish (see Fig. 3f), and the DAF order is not stabilized around this hole density.
Superconducting gap functions. The above mentioned electronhole asymmetry in the spin and the charge degrees of freedoms are the keys to understand the competition between the two types of SC. The main contribution of the gap function for the d xy -wave SC is given as where α(=1, 2) denotes a band index, and Δ α i is the pairing with the i-th neighbor dimers in the real space and treated as a variational parameter. We optimize the real space pairing up to 22nd neighbor dimers and find that the overall feature of the gap function is determined within the fourth neighbor, i.e., Δ α m for m ≤ 4. The term contaning Δ α 1 in Eq. (3) gives nodes in the horizontal (along k x -axis) and vertical (along k y -axis) directions. This is because the two diagonal pairings (orange and blue bars in Fig. 4a) have different sign, giving a d xy -type contribution. Therefore, this gap symmetry is referred to as d xy -wave 26 . Note that the terms corresponding to real space pairings parallel to horizontal (along x-axis) and vertical (along y-axis) directions vanish since they are along the nodal directions. This SC phase has been discussed in analogy with that of high-T c cuprates, ascribing the diagonal directions in κ-type structure to an approximate square lattice 10 . As in the case for high-T c cuprates, the sign of the gap function on the FS changes four times (see Fig. 4c).
On the other hand, the main contribution of gap function for the extended-s+d x 2 Ày 2 -wave SC is given as The first term contaning Δ α 1 in Eq. (4) (blue bars in Fig. 4b) does not change sign within the first Brilloin zone and becomes zero only along the zone boundary; this term gives an extended-slike contribution. Furthermore, Δ α 1 changes sign between different bands, namely, sgnΔ 1 1 ¼ ÀsgnΔ 2 1 . In this respect, the pairing symmetry can also be referred to as s ± , similar to that of ironbased SC 57,58 . The second and third terms containing Δ α 2 and Δ α 3 in Eq. (4) (orange and yellow bars in Fig. 4b, respectively) give nodes in the diagonal direction; these terms give a d  Fig. 1b). High symmetry points are Γ(0,0), M(π/R x , π/2R y ), X(π/R x ,0), and Y(0, π/2R y ). d-f Interdimer spin structure factor S dim (q) for (d) n hole = 2/3 with The gap symmetry is thus referred to as an extended-s + d x 2 Ày 2 -wave 26,28,32,59 . The sign changes of the gap function on the FS is shown in Fig. 4d. The competition between the two types of SC has been already discussed for the undoped case. First, in the effective 1/2-filled dimer model, many early studies inferred the d xy -wave type [12][13][14][16][17][18] . However, its stability over the AF phase has recently been doubted in several numerical studies including the VMC approach 20,22 . Kuroki et al. 26 were the first to point out the importance of treating the 3/4-filled model, namely, considering the intradimer charge degree of freedom. In the 3/ 4-filled Hubbard model, the two types of SC compete and either of them is favored depending on the parameters, especially the degree of dimerization 26,28,32 . It is only recently that the importance of the intersite Coulomb interaction V ij on SC was pointed out 27,29 . In fact, our previous VMC study for the undoped case found that, although both symmetries are enhanced by the intersite Coulomb interactions, the extendeds + d x 2 Ày 2 -wave is slightly more favored 29 .
As shown in Fig. 2f, our calculations find that this competition is released when the mobile carriers are doped into the system. For the hole-doped side, the AF spin fluctuation toward the DAF order is enhanced due to the shape of the FS and the Coulomb interactions, as already seen in Fig. 3, and consistent with the recent calculations based on the dimer model 46 . Similarly to the high-T c cuprates, the AF spin fluctuation mediated SC is then developed with the strong singlet correlation along diagonal bonds shown in Fig. 4a, resulting in the d xy -wave symmetry. On the other hand, the AF spin fluctuation is suppressed with electron doping as seen in Fig. 3f and the spin singlet correlation along horizontal and vertical bonds compete with that along diagonal bonds. The extended-s + d x 2 Ày 2 -wave is eventually favored since all bonds can contribute to the singlet pairing. The carrier doping deforms the shape of the FS and modifies the AF spin fluctuation, thus inducing the change of the symmetry of SC. Although the electron-hole asymmetry shown in Fig. 2a appears similar to that of high-T c cuprates at a glance, they are much different especially for the SC phase. In both cases, the van Hove singularity appears only in the hole-doped side, causing the electron-hole asymmetry. However, different physics are delicately involved in κ-(ET) 2 X, as we have shown so far, and even the change of the symmetry of SC occurs by carrier doping. This is due to the unique geometrical frustration inherent in the triangular-like lattice structure of κ-(ET) 2 X. Furthermore, this asymmetry is expected to be robust for κ-(ET) 2 X in general because the van Hove singularity is always located in the holedoped side for the realistic parameter set of these conductors.
We can show more directly how the spin and charge correlations are correlated to the stability of the SC phases. The interdimer charge structure factor is defined as where n dim l ¼ ðn 2lÀ1" þ n 2l" Þ þ ðn 2lÀ1# þ n 2l# Þ is the total charge density within l-th dimer formed by molecular sites 2l − 1 and 2l with the central position r l 56 . Figure 5a-d show S dim (q) and N dim (q) for n hole = 640/1152 = 0.556 and 544/1152 = 0.472, where the d xy -wave and the extended-s + d x 2 Ày 2 -wave SC are stabilized, respectively. For n hole = 0.556, S dim (q) peaks around (0, ±π/R y ) that are favorable for the d xy -wave SC, while for n hole = 0.472, S dim (q) shows frustrated spin structure that are favorable for the extended-s + d x 2 Ày 2 -wave SC. These are consistent with the mechanism of SC described above. On the other hand, N dim (q) peaks around (0, ±3π/2R y ) for both n hole = 0.556 and 0.472. This is because they are located near the instability of 3-fold CO-1 and 2, respectively, and the corresponding wave vectors are the same.
The superconducting condensation energy ΔE and spin/charge correlations are contrasting between the two SC phases. Figure 5e, f show the U=t b 1 dependence of N 3-fold , S peak , and ΔE. N 3-fold = N dim (0, ±3π/2R y ) and S peak = S dim (0, ±π/R y ) for n hole = 0.556. The peak position in S dim (q) changes along the upper and lower edges of the Brillouin zone for n hole = 0.472 and we take the maximum value for S peak . For n hole = 0.556, ΔE does enhance, but despite the more rapid increase of N 3-fold toward 3-fold CO-1 instability, it rather follows S peak . This is consistent with the usual AF spin fluctuation picture discussed in high-T c cuprates. On the other hand, for n hole = 0.472, ΔE is greatly enhanced following the rapid increase of N 3-fold , which suggests the close correlation between the extended-s + d x 2 Ày 2 -wave SC and the 3-fold CO-2 type charge fluctuation.

Discussion
Let us note that the intersite Coulomb interactions V ij are indispensable to the stability of SC, not only for the extended-s + d x 2 Ày 2 -wave SC phase but also for the metastable d xy -wave SC phase. As we have shown in the previous work 29 , no long-range ordered phases are stabilized in the absence of V ij for the undoped condition (or unphysically large U is necessary). This is also the case in the doped condition studied here. This indicates that in κ-(ET) 2 X, the charge degree of freedom is still active even with large U and the cooperation between U and V ij induces various phases such as the DAF, PCO, 3-fold COs, and SC. Especially, the extended-s + d x 2 Ày 2 -wave SC is enhanced toward both polar and 3-fold type CO instabilities. Recent experiment on photoinduced phase transition in κ-Br suggests that the polar charge oscillation is enhanced near the superconducting transition 60 , consistent with our picture that the SC is enhanced toward the CO instability. Experimentally, the symmetry of the SC for the undoped case is still controversial 38,[61][62][63][64][65][66][67] . This is consistent with our result suggesting that the competition between the two types of SC is most pronounced near the undoped n hole = 1/2. Slight modification of parameters may alter the stability of the two. Now in κ-Cl, both electron and hole doping are realized using the EDLT 46,47 . The filling-control MI transition and the emergence of SC have been confirmed with significant electron-hole doping asymmetry. Our result suggests that the different SC appears by carrier doping; the d xy -wave SC on the hole-doped side and the extended-s + d x 2 Ày 2 -wave SC on the electron-doped side. The former is due to the strong AF spin fluctuation, similar to high-T c cuprates, and the latter is the consequence of frustrated spin structure and enhanced charge correlation under the geometrical frustration characteristic of the κ-type ET compounds. Experiments for the chemically doped κ-(ET) 4 Hg 3−δ Y 8 suggest similarities with the high-T c cuprates, especially for the non-Fermiliquid behaviors above the SC transition temperature [42][43][44] . This is overall consistent with our results because this system is holedoped. However, the tight-binding parameters for this system are rather closer to the isotropic triangular-lattice like arrangement of dimers, where the AF phase is expected to be unfavored because of the stronger geometrical frustration. Indeed, the temperature dependence of the spin susceptibility is similar to that of κ-CN 45 , which does not show any magnetic long-range orders. It is difficult to fully explain the character of κ-(ET) 4 Hg 3−δ Y 8 within the present study and the analysis with appropriate tight-binding parameters is necessary for further understanding.
The nonmagnetic (candidate of gapless spin-liquid) insulator and neighboring SC observed in κ-CN are also an intriguing phenomena for the undoped case. Although the tight-binding parameters and the resulting electronic structure of κ-CN is different from the present study, the intradimer charge degree of freedom and the intersite Coulomb interactions should play crucial roles as pointed out previously 29 . Indeed, the anomalous dielectric response 36 and Raman spectroscopy 37 indicate that the intradimer charge degree of freedom is active in κ-CN. The stability of spin-liquid and SC phase in the 1/2-filled Hubbard model on the anisotropic triangular lattice, which is the approximate dimer model of κ-CN, is still controversial despite the long and extensive studies [12][13][14][15]17,[19][20][21][22][23][24][25] . The analysis for the 3/4-filled EHM employed in this study will be an alternative way to investigate this issue and it is left for future studies.

Methods
Details of the VMC method. Here, we show the details of the VMC method. Trial wave function is a Gutzwiller-Jastrow type, Ψ j i ¼ P J c P J s Φ j i. |Φ> is a one-body part constructed by diagonalizing the one-body Hamiltonian including the off-diagonal elements {D}, {M}, and {Δ} to treat long-range orders of charge, spin, and SC, respectively. The renormalized hopping integrals (t b 1 ;t b 2 ;t p ;t q ) are also included in |Φ> as variational parameters, wheret b 1 ¼ t b 1 is fixed as a unit. P J c ¼ exp ½À P i;j v c ij n i n j and P J s ¼ exp ½À P i;j v s ij s z i s z j are charge and spin Jastrow factors which control long-range charge and spin correlations, respectively. Here, s z i ¼ n i" À n i# , and we assume v c ij ¼ v c ðjr i À r j jÞ and v s ij ¼ v s ðjr i À r j jÞ, where r i is the position of molecular site i. The variational parameters in |Ψ> aret b 2 ,t p ,t q , {D}, {M}, fΔg, fv c ij g, and fv s ij g, and they are simultaneously optimized using the stochastic reconfiguration method 68 . The total number of molecular sites is 4 × L × L/ 2 = 2L 2 and varied from L = 12 to L = 24 with antiperiodic boundary conditions in both x and y directions of the primitive lattice vectors (see Fig. 1a).
The one-body part |Φ> for the PM, DAF, and PCO states can be obtained by diagonalizing the one-body Hamiltonian, with the hopping matrix elements given as T 32 ¼ À2t p e Àik y ðR y Àδ y Þ cos where c y kmσ (c kmσ ) is a creation (annihilation) operator of electron at molecule m (=1-4), as indicated in Fig. 1a, with momentum k and spin σ(=↑,↓), and s σ = 1 (−1) for σ = ↑(↓). M z m is a variational parameter which induces the staggered AF long-range order aligned to the z direction for molecule m; M z for the DAF state, and M z 1 ¼ M z 3 ≠M z 2 ¼ M z 4 for the PCO state. a y kασ (a kασ ) in Eq. (7) is a creation (annihilation) operator of quasiparticle in band α with momentum k and spin σ(=↑,↓), obtained by diagonalizing Eq. (6), andẼ α ðkÞ is the corresponding quasipaticle energy.
For the 3-fold CO-1 and 2 states, the unit cell is three times larger than the original one and contains 12 molecules, as shown in Fig. 2d, e. Therefore, |Φ> is obtained by diagonalizing a 12 × 12 matrix. The corresponding Weiss field is D mσ P kσ c y kmσ c kmσ , which induces charge disproportionation and spin ordering within the unit cell through the variational parameters D mσ for m = 1, 2, …, 12. Finally, |Φ> for SC is obtained by diagonalizing the BCS-type mean-field Hamiltonian, H BCS ¼ P k a y k1" ; a y k2" ; a y k3" ; a y k4" a Àk1# ; a Àk2# ; a Àk3# ; a Àk4# where Δ α denotes gap function for band α and ξ α ¼Ẽ α ðkÞ Àμ is a quasiparticle energy measured from the renormalized chemical potentialμ. Δ α ¼ Δ α ðfΔ α i gÞ is constructed from the real space pairing Δ α i up to the 22nd neighbor dimers (i = 1, 2, …, 22). Since the band 3 and 4 are located much below the Fermi energy (see Fig. 1b), their contribution to the pairing is negligible and therefore we set Δ 3 = Δ 4 = 0.

Data availability
The data that support the findings of this study are available from the corresponding author on reasonable request.

Code availability
The code that support the findings of this study are available from the corresponding author on reasonable request.