Two-dimensional ground-state mapping of a Mott-Hubbard system in a flexible field-effect device

See allHide authors and affiliations

Science Advances  10 May 2019:
Vol. 5, no. 5, eaav7282
DOI: 10.1126/sciadv.aav7282


A Mott insulator sometimes induces unconventional superconductivity in its neighbors when doped and/or pressurized. Because the phase diagram should be strongly related to the microscopic mechanism of the superconductivity, it is important to obtain the global phase diagram surrounding the Mott insulating state. However, the parameter available for controlling the ground state of most Mott insulating materials is one-dimensional owing to technical limitations. Here, we present a two-dimensional ground-state mapping for a Mott insulator using an organic field-effect device by simultaneously tuning the bandwidth and bandfilling. The observed phase diagram showed many unexpected features such as an abrupt first-order superconducting transition under electron doping, a recurrent insulating phase in the heavily electron-doped region, and a nearly constant superconducting transition temperature in a wide parameter range. These results are expected to contribute toward elucidating one of the standard solutions for the Mott-Hubbard model.


The electron correlation in solids, or the Mott physics, offers intriguing phenomena in materials such as unconventional superconductivity (SC). The key parameters for controlling the Mott physics are the electronic bandfilling and bandwidth. Variation of the former triggers the doping-induced high–transition temperature (high TC) SC in cuprates (1) and the field-induced SC in a twisted graphene bilayer (2), while variation of the latter triggers the pressure-induced SC in organics (3) and the chemical pressure–induced SC in fullerenes (4). However, the simultaneous control of these two parameters has been lacking so far, leaving a comprehensive phase diagram of correlated materials inaccessible. For example, high-TC cuprates can exhibit SC only in the bandfilling-controlled regime, and the bandwidth cannot be sufficiently controlled to induce SC from a nondoped Mott insulating state because of the hard crystal lattice. Recently, however, we have achieved both electric field–induced SC using a solid-gate field-effect transistor and strain-induced SC by substrate bending in similar organic Mott insulator devices (5, 6). This situation has motivated us to elucidate the two-dimensional ground-state map of the Mott-Hubbard model by combining these two technologies in a single organic device. Here, we report simultaneous control of the bandfilling and bandwidth at an organic Mott insulator interface, where the details of the superconducting transitions in the proximity of the Mott insulator in the two-dimensional ground-state map are revealed.

The two-dimensional organic Mott insulator κ-(BEDT-TTF)2Cu[N(CN)2]Cl (hereafter referred to as κ-Cl) comprises alternating layers of conducting BEDT-TTF+0.5 cations and insulating Cu[N(CN)2]Cl counteranions (7). The conducting BEDT-TTF molecules are strongly dimerized and can be modeled as a single-band Hubbard model on an anisotropic triangular lattice with t′/t = −0.44, where t is the nearest-neighbor (interdimer) hopping and t′ is the next-nearest-neighbor hopping (Fig. 1A) (8). Similarities between the κ-type BEDT-TTF salts and high-TC cuprates, such as the proximity between antiferromagnetism and SC, and the unconventional metallic phase in the normal state, have been pointed out (9). In contrast to the high-TC cuprates, where the SC is induced with carrier doping by chemical substitution, κ-Cl has been studied in terms of the pressure-induced superconducting transition because of its sensitivity to pressure and the difficulty of chemical doping. However, recently developed techniques for field-effect doping have enabled the reversible and finely tuned doping to κ-Cl (10, 11). In previous work (11), measurements of transport properties and calculations based on cluster perturbation theory indicate strong doping asymmetry, where major pseudogaps open near the van Hove critical points under substantial hole doping, while a more non-interacting-like Fermi surface appears under substantial electron doping. The doping asymmetry for the SC, which should be observed under pressure at lower temperatures, is currently of great interest because it will provide further insight into unconventional superconductivities in strongly correlated systems. Here, we report a two-dimensional mapping of the ground state for a Hubbard system based on an organic field-effect device. This is, to our knowledge, the first direct derivation of a ground-state phase diagram with fine control of the bandfilling and bandwidth in a Hubbard system in the solid state (Fig. 1B).

Fig. 1 Bandfilling and bandwidth control of κ-Cl in the same sample.

(A) Molecular arrangement of the BEDT-TTF layer in κ-Cl (top view). (B) Conceptual phase diagram based on the Hubbard model (37). The vertical axis denotes the strength of the electron correlation. κ-Cl is originally located near the tip of the insulating region and is shifted along both directions to investigate the superconducting region. (C) Schematic side view of the device structure. The doping concentration and effective pressure are controlled by an electric double-layer gating and bending of the substrate with a piezo nanopositioner, respectively. The resistivity is measured by the standard four-probe method. (D) Sheet resistivity versus temperature plots under hole doping at tensile strain S = 0.41%. The dashed line indicates the pair quantum resistance h/4e2. (E) Resistivity versus temperature plots under different tensile strains at gate voltage VG = 0 V.

To this end, we fabricated electric double-layer transistors (EDLTs) using thin single crystals of κ-Cl on flexible polyethylene terephthalate (PET) substrates (Fig. 1C and fig. S1). EDLT doping is particularly effective for materials with low carrier density because a small gate voltage can markedly vary the bandfilling in these materials. The carrier density of κ-Cl is much lower than that of high-TC cuprates (κ-Cl: ~1.8 × 1014 cm−2; YBa2Cu3O7−δ: ~6.7 × 1014 cm−2). According to the Hall effect and charge displacement current measurements, the carrier density in κ-Cl can be reversibly tuned up to approximately ±20% (Supplementary Text and fig. S2). In addition, the effective pressure of the sample can be simultaneously tuned by applying a bending strain to the substrate (6). Therefore, this combination of the EDLT doping and bending-strain techniques allows us to investigate the doping (δ)–correlation (U/t)–temperature (T) phase diagram for the same sample (U: on-dimer coulomb repulsion).


First, we study the resistivity without gating. Figure 1E shows the temperature dependence of the resistivity ρ for various values of the tensile strain S at zero gate voltage. As can be seen, the resistivity of κ-Cl is highly sensitive to the strain because it is in close proximity to the pressure-induced first-order Mott insulator/superconductor transition (12, 13). κ-Cl, which is originally a Mott insulator at ambient pressure, is a superconductor on a PET substrate owing to the compressive strain from the substrate induced by its thermal contraction at low temperatures. Upon applying the tensile strain S, κ-Cl exhibits a strain-induced superconductor-to-insulator transition at the lowest temperature. At higher temperatures, the resistivity shows a nonmonotonic dependence, changing from semiconducting to metallic with decreasing temperature, followed by resistivity jumps to antiferromagnetic Mott insulating states at lower temperatures. This peculiar dependence on the temperature indicates that the strained κ-Cl is in the proximity of the superconductor/antiferromagnetic Mott insulator transition at half filling, similarly to the bulk sample under hydrostatic pressure (1315).

Next, we study the resistivity under gating. For S ≤ 0.44%, ambipolar SC is observed under gating (Figs. 1D and 2). Upon applying negative gate voltages, the resistivity monotonically decreases by orders of magnitude and SC emerges for VG ≤ −0.3 V at S = 0.41% (Fig. 2C). The critical temperature TC of approximately 12 K does not vary notably down to VG = −0.5 V.

Fig. 2 Electron-hole asymmetric phase diagram of κ-Cl.

(A to G) Contour plots of the sheet resistivity ρ under tensile strains S = 0.35% (A), 0.39% (B), 0.41% (C), 0.44% (D), 0.46% (E), 0.50% (F), and 0.55% (G) as a function of temperature and gate voltage. (H) Contour plots of the sheet resistivity ρ at 5.5 K as a function of gate voltage and tensile strain. Black dots in all figures indicate the data points where the sheet resistivity was measured. The doping concentration estimated from the average density of charge accumulated in the charge displacement current measurement (fig. S2) is shown for reference on the upper horizontal axis in (H). AFI, antiferromagnetic insulator; h-SC, p-type SC; e-SC, n-type SC. In the region below the white dashed line at S = 0.35%, the surface resistivity under doping cannot be measured because the nondoped bulk is superconducting below 12 K.

On the other hand, the effect of VG is not monotonic for electron doping. The resistivity abruptly drops, and a superconducting state with TC of approximately 12 K again emerges with low electron doping (+0.14 V ≤ VG ≤ +0.22 V at S = 0.41%). However, the resistivity increases again, and the SC disappears after further electron doping. Although the hole- and electron-doped superconducting states show similar values of TC, the doping concentration where SC appears is significantly doping asymmetric, as shown in Fig. 2 (B to D). Because the charge neutrality point (highest resistivity point) corresponding to half filling is approximately +0.05 V in this sample, the electron-doped SC is observed only in a narrow region of VG near half filling. These superconducting states are suppressed by magnetic fields, and their normal states are insulating (fig. S3), as in the pressure-induced superconducting state at half filling (16). Curiously, the normal-state resistivity takes a local minimum against VG at the optimum doping on the electron-doped side, as shown in the upper part of Fig. 2 (B to D). The “dip” of the resistivity is observed at least up to 80 K, which is several times higher than TC (fig. S4). The superconducting regions in the VG-T plots become narrower and eventually disappear with increasing tensile strain (Fig. 2, F and G). Both the p- and n-type SC fade at almost the same value of strain.

Summarizing the tensile strain dependence of ρ in the VG-T plots (Fig. 2, A to G), we obtain the VG-S plots at 5.5 K (Fig. 2H), which should reflect the ground-state bandfilling-bandwidth phase diagram. (However, note that the sample conductance is the sum of the doped surface conductance and the nondoped bulk conductance. That is, the sample resistivity deviates from the doped surface resistivity unless the bulk conductance is much less than the surface conductance. In particular, the resistivity data at S = 0.35% below 12 K do not reflect the doped states. Once the nondoped bulk becomes superconducting, it cannot be determined whether the doped surface is superconducting.)

While the hole-doped superconducting phase lies slightly away from the antiferromagnetic insulating phase, the electron-doped superconducting phase is located in a very narrow region next to the antiferromagnetic insulating phase. The resistive transition between the half-filled insulating state and the electron-doped superconducting state is very abrupt (ρ falls by up to eight orders of magnitude within a gate voltage range of ~0.1 V). At the boundary, discontinuous fluctuations in resistivity are also observed (fig. S5). These behaviors are reminiscent of the pressure-induced phase separation/percolative SC at half filling (1315). Therefore, it is reasonable to consider that the first-order transition line at least extends from half filling to the electron-doped superconducting phase boundary in the diagram.

To summarize the experimental results, (i) the superconducting region is narrower and closer to half filling under electron doping than under hole doping, (ii) the electron doping–driven superconducting transition is very abrupt and discontinuous (first-order–like), (iii) the superconducting regions appear to be connected to each other in the VG-S plots within the data resolution, and (iv) the doping asymmetry for the resistivity diminishes with decreasing bandwidth.

In contrast to other solid-state materials such as the high-TC cuprates and doped fullerenes, the bandfilling and bandwidth are simultaneously controlled in a single sample, and thus, the sample dependence is less significant here. Furthermore, one simple molecular orbital (highest occupied molecular orbital) governs the electronic properties under both electron and hole doping in κ-Cl. Accordingly, the observed doping asymmetry is attributed to the intrinsic electron-hole asymmetry of interacting fermions on the anisotropic triangular lattice.

For a better understanding of the particle-hole asymmetry, we used the variational cluster approximation (VCA) (17) to theoretically consider antiferromagnetic and dx2y2 superconducting orders in the Hubbard model on the anisotropic triangular lattice defined by the following HamiltonianH^=ij,σtij(c^iσc^jσ+H.c.)+Uin^in^iμiσn^iσ(1)where c^iσ creates an electron site i with spin σ(,),n̂iσ=ĉiσĉjσ,tij is the transfer integral between neighboring sites i and j (indicated as t and t′ in Fig. 1A), U is the on-site coulomb repulsion, and μ is the chemical potential. In the following, we set t′/t = −0.44 (8) and used a cluster of size L = 4 × 3. We performed scans at 0 K, varying μ for several values of U/t. The strain dependence of the intradimer charge degree of freedom is not considered because the effective coulomb repulsion on the dimer (U) is not sensitive to changes in the intradimer transfer integral (5, 18): Because U2|tid|4tid2/Ubare, where tid is the intradimer transfer integral (~0.25 eV) and Ubare is the on-site coulomb repulsion for a single BEDT-TTF molecule (~1.0 eV), small changes in tid cancel out. Here, we used the one-band model to obtain an approximate outline of the phase diagram. However, note that it is not perfect for describing the electronic properties of κ-type BEDT-TTF salts. Different gap symmetry (extended s- + d-wave symmetry) has been predicted using the four-band models (1921) and observed by scanning tunneling microscopy in a superconducting salt (22).

Figure 3A shows the antiferromagnetic and superconducting order parameters as a function of doping concentration δ = n − 1, where n is the electron density and δ = 0 corresponds to half filling. Similarly to the experiment, the antiferromagnetic order breaks down, and SC emerges at higher δ. Even with a small amount of electron doping (δ < 0.1), the transition takes place when U/t is small. In such a case, the SC appears more abruptly than for hole doping. However, this abrupt breakdown of the antiferromagnetic insulating phase disappears when U/t is increased to more than 5.0.

Fig. 3 VCA calculations.

(A) Antiferromagnetic and dx2y2 superconducting order parameters, M and D, respectively, versus doping concentration δ for several values of U/t. M and D for the metastable and unstable solutions (empty symbols) are also shown at U/t = 4 and 4.5 under electron doping (corresponding to positive δ). (B) Doping concentration δ versus chemical potential μ relative to that at half filling (μhalf) for several values of U/t. The results for the metastable and unstable solutions at U/t = 4 and 4.5 are indicated by dashed lines, while the results obtained by the Maxwell construction are denoted by solid vertical lines. This implies the presence of phase separation and a first-order phase transition. It is noteworthy that there is a steep (nearly vertical) increase in δ with increasing μ for larger values of U/t under electron doping, suggesting a strong tendency toward phase separation. The values of μhalf are μhalf = 1.3725t for U/t = 3.5, μhalf = 1.8375t for U/t = 4, and μhalf = U/2 for U/t 4.5. (C) Chemical potential μ versus doping concentration δ for U/t = 4 (see fig. S8 for more details). δ1 and δ2 are the doping concentrations of the two extreme states in the phase separation. All results in (A) to (C) are calculated using the VCA for the single-band Hubbard model on an anisotropic triangular lattice (t′/t = −0.44) with a 4 × 3 cluster. (D) Noninteracting tight-binding band structure and density of states (DOS) for t′/t = −0.44 with t = 65 meV. Here, Γ = (0,0), Z = (0,π/c), M = (π/a,π/c), and X = (π/a,0), with a and c being the lengths of the primitive translation vectors indicated in Fig. 1A. The Fermi level for half filling is set to zero and denoted by dashed lines. (E) Single-particle spectral functions and DOS of κ-Cl at half filling in an antiferromagnetic state at zero temperature, calculated by VCA. The Fermi level is denoted by a dashed line at zero energy. The flat features seen away from the Fermi level indicate incoherent continuous spectra due to the electron correlation. The reason why they appear rather discretized is because of the discrete many-body energy levels in the VCA calculation, for which a finite-size cluster is used to obtain the single-particle excitation energies.

In addition, we found a nonmonotonic dependence of the chemical potential on δ under electron doping (Fig. 3B). Accordingly, Fig. 3A includes metastable and unstable solutions on the electron-doped side for small U/t (indicated as empty symbols). As shown in Fig. 3C, the Maxwell construction reveals phase separation between two phases with different doping concentrations δ1 and δ2, where the volume fraction of each phase is proportional to δ − δ1 or δ2 − δ for the average doping concentration δ, assuming that δ1 ≤ δ ≤ δ2. If one of the phases is superconducting, percolation SC can appear when δ exceeds the percolation limit. Although the results depend quantitatively on the clusters used, we found that the tendency toward phase separation is robust under electron doping (see section S2).

According to the resistivity behavior without a gate voltage (fig. S6), the change in U/t from S = 0.35 to 0.55% is comparable to or slightly larger than that from bulk crystals of κ-(BEDT-TTF)2Cu[N(CN)2]Br (κ-Br, U/t = 5.1, superconducting) to κ-Cl (U/t = 5.5, Mott insulating) (8). That is, the variation in the experiments (~8% change in t for 0.2% increase in strain) is considered to be smaller than that in the calculations shown in Fig. 3A. In addition, the calculations indicate the Mott transition at half filling between U/t = 3.5 and 4.0, while it occurs between κ-Br (U/t = 5.1) and κ-Cl (U/t = 5.5) in the real material system. Despite these discrepancies, the characteristic features in the proximity of the transitions such as the doping asymmetry and the phase separation on the electron-doped side are qualitatively reproduced by our calculations.

The tendency of the phase separation with the sign of doping is opposite to that in the Hubbard model on a square lattice with the model parameters relevant for cuprates, where the phase separation is more likely to occur under hole doping (23). We attribute the phase separation tendency to the high density of states around the bottom of the upper Hubbard band split off from the van Hove singularity by the electron correlation (Fig. 3E). The high density of states accumulated along the Z-M line is characteristic of the anisotropic triangular lattice owing to the fact that the direction of the hopping t′ is orthogonal to the Z-M line, along which the dispersion remains flat. The calculations imply that our experiment was able to capture such a highly nontrivial correlated band structure on the unoccupied side as an abrupt change in the transport properties under electron doping.


The doping and correlation dependence of TC in the Hubbard model on an anisotropic triangular lattice has also been calculated using cluster dynamical mean field theory (24). It predicts hole-doped superconducting states in a wide doping region from extremely low to moderate hole doping (from 1 to 10%), while no SC appears under low electron doping of less than 10%. TC tends to be enhanced under hole doping and reduced under electron doping, compared with the value at half filling. However, in our experimental results shown in Fig. 2, there are many discrepancies among these theories as well as our previous expectations (5). An abrupt first-order superconducting transition was expected, but it is only observed in the electron-doped region. The recurrent insulating phase in the heavily electron-doped region is also a new observation. The values of TC are unexpectedly independent of both doping and strain. In addition, the electron-doped SC seems to persist at higher U/t than the hole-doped SC. Note that the optimal doping concentration does not change significantly when the strain is varied. These are all new findings observed by the two-dimensional scanning of the strain and carrier density at the κ-Cl surface, although their origins require further discussion.

The disappearance of SC and the appearance of an insulating phase for large electron doping (VG ≥ +0.25 V) cannot be explained by the calculations, for example. To understand this, we should probably be aware of the simplicity of the model considered here, which ignores elements specific to this material, such as the intradimer charge degree of freedom and the intersite coulomb repulsion (20). One possibility is that a magnetic or charge-ordered state emerges at specific doping levels (for example, ~12.5%), as in the case of the stripe order in the cuprates (2528). Recent theoretical calculations for the Hubbard model on a square lattice suggest that a large part of the macroscopic phase separation region can be replaced by more microscopically inhomogeneous stripe states (29, 30). The presence of a magnetic or charge-ordered state with high resistivity at 12.5% doping may explain the dip in resistivity in the normal state (fig. S4).

The correlation strength and bandfilling are the most fundamental parameters that determine unconventional SC in correlated electron systems. To construct an experimental correlation-bandfilling phase diagram, however, it has been necessary to combine data from different materials, resulting in unavoidable material and sample dependences. Obtaining the phase diagram in a single sample is a milestone for understanding the pristine mechanism of unconventional SC. This is achieved here by virtue of the high tunability of the organic Mott insulator. In addition, despite the use of a real material, the high controllability and cleanness of our device are in common with those of quantum simulators (31) such as cold atoms in an optical lattice. We hope that these results also serve as a useful reference for quantum simulations of the Fermi-Hubbard model beyond the field of materials science.


Sample preparation and transport measurement

The source, drain, and gate electrodes (18-nm-thick Au) were patterned on a PET substrate (Teflex FT7, Teijin DuPont Films Japan Limited) using photolithography. A thin (~100 nm) single crystal of κ-Cl was electrochemically synthesized by oxidizing BEDT-TTF (20 mg) dissolved in 50 ml of 1,1,2-trichloroethane [10% (v/v) ethanol] in the presence of TPP[N(CN)2] [tetraphenylphosphonium (TPP), 200 mg], CuCl (60 mg), and TPP-Cl (100 mg). After applying a current of 8 μA for 20 hours, the thin crystal was transferred into 2-propanol with a pipette and guided onto the top of the substrate. A diagonal of the rhombic crystal, which is usually parallel to the crystallographic a or c axis, was aligned parallel to the direction of the current and strain (although it is unknown which axis was the a or c axis). After the substrate was removed from the 2-propanol and dried, the κ-Cl crystal was shaped into a Hall bar using a pulsed laser beam with a wavelength of 532 nm. The typical dimensions of the Hall bar sample were approximately 15 μm (width) by 40 μm (length) by 100 nm (thickness), as shown in fig. S1. Figure S1C shows an atomic force microscopy image of the surface of a typical κ-Cl crystal laminated on the substrate. The roughness of the surface was suppressed to less than 1.5 nm over the micrometer scale (fig. S1D), which corresponds to the thickness of one set of BEDT-TTF and anion layers. As a gate electrolyte, the ionic liquid 1-ethyl-3-methylimidazolium 2-(2-methoxyethoxy) ethylsulfate was added dropwise to the sample and an Au side gate electrode. The EDLT device was completed by mounting a 1.2-μm-thick polyethylene naphthalate (PEN) film on the ionic liquid droplet. Thinning of the gate electrolyte using the PEN film reduced the thermal stress at low temperatures.

The transport measurements were performed using a Physical Property Measurement System (Quantum Design). The four-terminal resistance was measured with a dc of 1 μA using a dc source (Yokogawa 7651, Yokogawa) and a nanovoltmeter (Agilent 34420A, Agilent Technologies). The maximum applied voltage was limited to 1 V (that is, the measurement was voltage biased in the high-resistance states). The applied current was monitored with a current amplifier (SR570, Stanford Research Systems). The resistivity was estimated from the resistance and sample dimensions. In Fig. 1E, the resistivity ρ (Ωcm) was estimated as ρ = R × W × D/L, where R, W, D, and L denote resistance, width, thickness, and length of the sample, respectively. In the other figures, the sheet resistivity ρ = R × W/L is shown because only the sample surface was doped and the three-dimensional resistivity cannot be accurately estimated. According to the Hall measurements (11), the effect of doping is confined within one molecular layer under sufficient doping.

The four-terminal resistivity was measured during temperature cycles between 220 and 5 K. The temperature was varied at rates of 2 K/min (T > 20 K) and 0.3 K/min (T < 20 K). The gate voltage VG was swept from +0.5 to −0.5 V with a step of 0.05 V at 220 K, but in the low–electron doping regime from +0.22 to 0.08 V, it was more finely tuned with a step of 0.02 V. After varying the gate voltage, we waited 1 min for the stabilization of the gate bias before cooling the sample.

Because the samples are top-gate and bottom-contact transistors, the current flows through the ungated regions between the electrodes and the doped surface, resulting in non-negligible contact resistance. Figure S1 (E and F) shows the four- and two-terminal resistances without and with a gate voltage, respectively. One can see that the temperature dependence of the contact resistance is more moderate than that of the four-terminal resistance. Therefore, the ratio of contact resistance to sample resistance is large when the sample is metallic or superconducting. However, the contact resistance is typically up to the order of 10 kilohms and does not hinder the four-terminal resistivity measurement.

After the temperature cycles at different gate voltages, a tunable tensile strain was mechanically applied with a nanopositioner (ANPz51, attocube) from the back side of the substrate at 220 K. The basic techniques and apparatus for the strain measurements were the same as those in our previous report (6). The tensile strain was estimated from S = 4tx/(l2 + 4x2), where t and l are the thickness and length of the substrate, respectively, and x is the displacement of the nanopositioner. We assumed that (i) the bent substrate (and κ-Cl crystal) is an arc of a circle and (ii) the angle between the ends of the substrate is small (small-angle approximation). The tensile strain was applied in the following order: 0.39, 0.41, 0.44, 0.46, 0.50, 0.55, and 0.35%.

VCA calculations

1. Model Hamiltonian. To study the ground-state phase diagram of the organic Mott insulator κ-Cl from the theoretical point of view, we considered the single-band Hubbard model on the anisotropic triangular lattice (8, 3235). The Hamiltonian of the model is given by Eq. 1. The transfer integral between the different (same) dimers of BEDT-TTF molecules is given by tij = t (t′). We set t′/t = −0.44, which is relevant for κ-Cl (8), and varied U and μ to control the strength of the electron correlation and the carrier concentration, respectively. In experiment, U/t and μ can be controlled by applying the strain S and the gate voltage VG, respectively. That is, the tensile strain S increases U/t, and the positive (negative) gate voltage increases (decreases) μ.

2. Variational cluster approximation. We applied the VCA (36) to the single-band Hubbard model in Eq. 1. The VCA is a many-body variational method based on the self-energy functional theory (3739) and allows us to investigate the possible spontaneous symmetry breaking including antiferromagnetism (AFM) and SC by taking into account the short-range electron correlations through the exact self-energy of a so-called reference system, which will be introduced later. The VCA has been applied for the single-band Hubbard model on the anisotropic triangular lattice at half filling (4042) and also on the square lattice with finite dopings relevant for cuprates (23, 43) to study possible magnetism and SC. Moreover, the VCA has been able to capture the first-order phase transtions of strongly correlated systems, as well as the free-energy balance and the Maxwell construction consistently (4446).

Let us first review the variational principle for the grand potential on which the VCA was constructed (38, 47, 48). There exists a functional Ω[Σ] of the self-energy Σ that gives the grand potential Ω as Ω = Ω[Σ*] such thatδΩ[Σ]δΣ|Σ=Σ*=0(2)whereΩ[Σ]=F[Σ]1βTr ln(G01+Σ)(3)G0 is the noninteracting single-particle Green’s function, F[Σ] is the Legendre transform of the Luttinger-Ward functional Φ[G] (45), i.e., F[Σ] = Φ[G] − β−1Tr[GΣ] with Σ = βδΦ[G]/δG, G is the interacting single-particle Green’s function, β = 1/T is the inverse temperature, and Tr is the functional trace that runs over all (both spatial and temporal, regardless of discrete or continuous) variables of the functions in Tr[· · ·]. The stationary condition Eq. 2 ensures that Σ* is the physical self-energy in the sense that Σ* satisfies the Dyson equation. The variational principle for the grand potential has been derived by the diagrammatic-expansion technique for many-body Green’s function (47) and later by the functional-integral technique for the grand-partition function in a nonperturbative way (48). A remarkable property of the Luttinger-Ward functional is that its functional form F depends only on the interaction term of the Hamiltonian (48).

Although the variational principle is rigorous and nonperturbative, a difficulty in applying it for practical calculations is that the explicit form of the Luttinger-Ward functional F[Σ] is unknown. The key approximation of the VCA is to restrict the space of the trial self-energy Σ to that of the exact self-energy Σ′ of a reference system described by the Hamiltonian Ĥ′ whose interaction term must remain unchanged, but the single-particle terms can be modified from Ĥ for Ĥ′ to be solvable exactly. The invariance of the Luttinger-Ward functional F enables us to write the grand-potential functional of reference system as Φ[Σ]=F[Σ]β1Tr ln(G01+Σ), where G0 is the noninteracting Green’s function of the reference system. By restricting the trial self-energy to Σ′ and eliminating F[Σ] from Eq. 3, we obtainedΩ[Σ]=Ω[Σ]1βTr ln(IVG[Σ])(4)where V=G01G01 represents the difference of the single-particle terms between the original and reference systems and G[Σ]=(G01Σ)1 is the exact interacting single-particle Green’s function of the reference system. Note that the right-hand side of Eq. 4 is computable as long as the reference system Ĥ′ is solvable.

The remaining arbitrariness of the single-particle terms in the reference system Ĥ′ allows us to optimize the trial self-energy Σ′ through varying single-particle fields λ, which parameterize the single-particle terms, as variational parameters, i.e., Σ′ = Σ′(λ), so as to satisfy the stationary condition in Eq. 2. Simply denoting Ω(λ)=Ω[Σ(λ)], the optimization scheme of the VCA amounts to find the extremaΩ(λ)λ|λ=λ*=0(5)in the space of the variational parameters λ, where λ* denotes the set of optimal variational parameters.

3. Reference system. Now, we specify our reference system Ĥ′ used for the VCA calculations. We assumed that our reference system is composed of a collection of identical, disconnected, and finite-size clusters, each of which is described by a Hamiltonian Ĥc(Ric), i.e., H^=ic=1NH^c(Ric), where Ric is the icth position of the cluster and N is the number of the clusters in the reference system. Because the clusters are identical, i.e., H^c(Ric)=H^c, we refer to it as Ĥc.

To study the SC and AFM phases under the carrier doping, we considered the following Hamiltonian of the clusterH^c=H^H+H^ϵ+H^h+H^Δ(6)whereH^ϵ=ϵi(n^i+n^i)(7)H^h=hi(n^in^i)eiQri(8)H^Δ=Δi,j(ηijc^ic^j+H.c.)(9)and ĤH is the Hubbard Hamiltonian given in the right-hand side of Eq. 1 but defined in the cluster. Here, Q = (π, π) and ri is the position of the ith site. ϵ is the variational parameter for the on-site potential that has to be optimized to calculate the particle density satisfying the thermodynamic consistency (23). Δ and h are the variational parameters that are introduced to detect the d-wave SC and AFM states by explicitly breaking the corresponding U(1) and SU(2) symmetries in the reference system, respectively. The spatial dependence of the form factor ηij for the SC considered here is given asηij={+1if rirj=±e1,1if rirj=±e2,0 otherwise,(10)where e1 and e2 are the primitive lattice vectors on the anisotropic triangular lattice bridging the two sites connected by the nearest-neighbor hopping integral t (not by t′). Note that the form factor ηij corresponds to that of the dx2y2–wave SC if it is considered on the square lattice spanned by the primitive lattice vectors e1 and e2 (assuming that e1 and e2 point to the x and y directions, respectively). Therefore, we refer to the SC state found in the VCA as dx2y2–SC.

4. Particle-hole transformation. Because the pairing term ĤΔ in Eq. 9 explicitly breaks the U(1) symmetry, the particle number is not conserved as long as Δ is finite. This is inconvenient if the computer program of the exact diagonalization method is implemented with the number of particles fixed. However, because the z component of the total spin is conserved, it is still possible to work in the basis set for the fixed number of particles after the particle-hole transformation{c^i:=c^id^i:=c^i(11)

In terms of the newly defined spinless fermion operators ĉi and d̂i, the single-particle terms in Eqs. 7 to 9 becomeH^ϵ=ϵi(n^icn^id+1)(12)H^h=hi(n^ic+n^id1)eiQri(13)H^Δ=Δi,j(ηijd^ic^j+H.c.)(14)where n̂ic=ĉiĉi and n̂id=d̂id̂i. Note that ĤΔ becomes the hybridization between the c and d orbitals. By applying the particle-hole transformation also for ĤH, the Hamiltonian of the cluster is now given asH^c=ijtij(c^ic^j+H.c.)+ijtij(d^id^j+H.c.)+(Uμ+ϵ+h)in^ic+(μϵh)in^id(15)Uin^icn^id+Δij(ηijd^ic^j+H.c.)μL+ϵLhieiQriwhere L is the number of sites in a cluster. It is now apparent that the total number of the particles is conserved as the operator commutes with the Hamiltonian Ĥc, which is inherited from the conservation of the z component of the total spin before the particle-hole transformation.

Here, several remarks on the nonoperator terms (i.e., constant terms) in Eq. 15 are in order. Unlike the single-particle operator terms multiplied by the variational parameters, the nonoperator terms in the last line of Eq. 15 cannot be “subtracted” as perturbation in the V matrix in Eq. 4 because the V matrix is a representation of the single-particle operators (34). Instead, those numbers have to be subtracted either from Ĥc itself or from the eigenspectrum of the cluster Hamiltonian Ĥc (for example, see Eq. 16). Note also that the term −μL in the third line of Eq. 15 has to be kept because μ is the model parameter that has the definite physical meaning. Last, we comment on the number of the basis states in the cluster. Because we considered the ground state of the cluster at zero temperature, we can focus on the sector for which the z component of the total spin is zero, i.e., i=1L(n^in^i)|x=0|x, where |x〉 is an arbitrary basis state for the corresponding sector. After the particle-hole transformation of Eq. 11, this sector corresponds to that with the total number of particles being L, i.e., i=1L(n^ic+n^id)|x=L|x. Because the cluster consists of 2L orbitals with L particles, the number of configurations in this sector is given by the binomial coefficient (2LL). In the main text, we considered the cluster of the 4 × 3 sites (L = 12), for which the number of the configuration is (2412)=24!/(12!)2=2704156. Considering that the VCA has to calculate the ground state and the single-particle Green’s function of the cluster repeatedly until the stationary condition Eq. 5 is satisfied and also that the model parameters U/t and μ/t are swept finely in a wide range to find not only stable but also metastable and unstable solutions that can identify the regions of the phase separation or the first-order transitions, L = 12 is reasonably large among the clusters feasibly handled by the exact diagonalization method.

5. Grand-potential functional and order parameters. Having specified the reference system, we can further substantiate the grand-potential functional. Because the system is considered at equilibrium and the clusters in the reference system are periodically aligned, the summand of the functional trace Tr[· · ·] in Eq. 4 is diagonal with respect to the Matsubara frequencies and the wave vectors of the superlattice on which the clusters are lined up. The grand-potential functional per site can hence be written asΩ=Ω1NLβνkln det[IV(k)G(iων)](16)whereΩ=1Lβlnseβ(EsE)(17)

Es is the eigenvalue of Ĥc for sth eigenstate |Ψs⟩ with the ground state |Ψ0⟩ (note that Ĥc includes the chemical potential term in our convention), E=ϵLhieiQri is the nonoperator term in Ĥc as annotated after Eq. 15, ων = (2ν + 1)π/β is the fermionic Matsubara frequency with ν integer, and k is the wave vector defined in the Brillouin zone of the superlattice with the number of wave vectors k being N. Here, we redefined Ω and Ω′ as the corresponding quantities Ω[Σ]′ and Ω′[Σ′] per site. G′(z) represents the exact single-particle Green’s function of the cluster, which is given asGijcd(z)=seβ(ΩEs)(Gij,s+,cd(z)+Gij,s,cd(z))(18)whereGij,s+,cd(z)=Ψs|c^i[z(H^cEs)]1d^j|Ψs(19)Gij,s,cd(z)=Ψs|d^j[z+(H^cEs)]1c^i|Ψs(20)and similarly for Gijcc(z), Gijdc(z), and Gijdd(z). The block Lanczos method is adapted for the calculation of G′(z). V represents the hopping process between the clusters and the subtraction of single-particle terms added to the reference system (i.e., Ĥϵ, Ĥh, and ĤΔ after the particle-hole transformation without the nonoperator term).

For a given set of the model parameters U/t and μ/t, all the variational parameters λ = (ϵ, h, Δ) are optimized simultaneously to satisfy the stationary condition(Ω(λ)ϵ,Ω(λ)h,Ω(λ)Δ)|λ=λ*=(0,0,0)(21)where λ* = (ϵ*, h*, Δ*) are the optimal variational parameters. The Newton-Raphson method was used for the optimization. The first and the second derivatives of Ω with respect to the variational parameters necessary for the gradient and the Hessian matrix are calculated by the central finite-difference method with the error in the second order of the step size. The step size for the finite difference was chosen adaptively by monitoring the curvature of Ω in the variational parameter space according to the scheme proposed in (45).

After the optimal variational parameters λ* were found, the expectation value of a single-particle operator could be evaluated from the VCA Green’s functionG(k,iων)=G(iων)[IV(k)G(iων)]1|λ=λ*(22)

For example, the doping concentration δ, the staggered magnetization M, and the dx2y2–SC order parameter D were calculated, respectively, asδ=1NLβνkTr[G(k,iων)](23)M=12NLβνkTr[mG(k,iων)](24)D=12NBβνkTr[ηG(k,iων)](25)where B is the number of pairs of sites connected by the dx2y2–wave pairing field in the cluster, and m and η are real symmetric 2L × 2L matrices whose matrix elements are either −1, 0, or +1 to represent the staggered magnetization and the dx2y2 pairing, respectively. Note that G(k,iων) was assumed to be in the Nambu spinor representation, and thus, the trace of G(k,iων) is not the electron density but the doping concentration.

In the zero-temperature limit, the sum Σs over the eigenstates of Ĥc in Eqs. 17 and 18 was taken only for the ground state of the cluster (s = 0), and the sum Σν over the Matsubara frequencies in Eqs. 18 and 23 to 25 was converted to the integral over the continuous frequency along the imaginary axis with a proper regularization for the integrand (49). The number N of wave vectors k in the Brillouin zone of the superlattice was chosen adaptively by monitoring the convergence of the sum Σk over k with respect to N for each frequency. In general, the larger N is required for frequencies closer to the Fermi level to reach the convergence.


Supplementary material for this article is available at

Section S1. Estimation of injected charge density

Section S2. Phase separation tendency of the VCA calculations

Fig. S1. Optical images, surface profile, and contact resistance properties of the flexible EDLT based on κ-Cl.

Fig. S2. Estimation of injected charge density.

Fig. S3. Suppression of SC by applying a magnetic field.

Fig. S4. Resistivity dip under small electron doping.

Fig. S5. Resistance versus temperature plots under low electron doping at tensile strain S = 0.41%.

Fig. S6. Comparison of temperature dependences of the resistivity between our device and bulk crystals.

Fig. S7. Cluster size dependence of the phase separation tendency.

Fig. S8. Swallowtail-shaped grand potential and Maxwell construction under electron doping with moderate correlation.

Reference (50)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to acknowledge Teijin DuPont Films Japan Limited for providing the PET films. Funding: This work was supported by MEXT and JSPS KAKENHI (grant nos. JP16H06346, JP15K17714, JP26102012, JP25000003, and JP16H04140), JST ERATO, MEXT Nanotechnology Platform Program (Molecule and Material Synthesis), and MEXT HPCI Strategic Programs for Innovative Research (SPIRE) (grant nos. hp140128 and hp150112). We are also grateful for allocation of computational time of the HOKUSAI GreatWave and HOKUSAI BigWaterfall supercomputer at RIKEN Advanced Center for Computing and Communication (ACCC) and the K computer at RIKEN Center for Computational Science (R-CCS). K.S. acknowledges support from the Overseas Research Fellowship Program of the Japan Society for the Promotion of Science. Author contributions: Y.K. and K.S. contributed equally to this work. Y.K. performed the planning, sample fabrication, cryogenic transport measurements, and data analyses. K.S. and S.Y. performed all the VCA calculations and data analyses. J.P. and T.T. developed the techniques for ion liquid and improved the device performance. S.T. performed the displacement current measurements and the data analyses. Y.K., K.S. and H.M.Y. wrote the manuscript. T.T., S.Y., H.M.Y., and R.K. supervised the investigation. All authors commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article