Optimal charge inhomogeneity for the d+id-wave superconductivity in the intercalated graphite CaC₆

Abstract

The coexistence of superconductivity and charge inhomogeneity was observed in many cuprate superconductors. The relationship between those two is still controversial. Similarly, in the graphene sheets of the intercalated graphitic superconductor CaC₆, the charge inhomogeneity was also observed. We simulate such a system by constructing the Hubbard model on the honeycomb lattice with charge inhomogeneity imposed by force. Utilizing the finite-temperature determinant quantum Monte Carlo algorithm, we examine the relationship between the superconducting pairing and the charge inhomogeneity. An optimal charge inhomogeneity for the d+id-wave pairing is found. While for other artificial charge inhomogeneities, the d+id-wave pairing is monotonically suppressed. The possible π-phase shift induced by charge inhomogeneity is also examined.

Graphical abstract

Keywords

intercalated graphitic superconductor / charge inhomogeneity / Hubbard model / d+id-wave pairing

1 Introduction

Inhomogeneity plays an important role in the strongly-correlated systems [1]. One of the famous phenomena is the formation of stripes in cuprate superconductors, i.e., the distribution of charge inhomogeneity with a period of around $4$ lattice spacings [2–7]. Despite intensive studies, the relationship between superconductivity (and its $d$-wave pairing mechanism) and charge inhomogeneity is still controversial, both experimentally [8–11] and theoretically [12–19]. In particular, the interrelation is argued to be multifacetted, i.e., the superconducting state and charge inhomogeneity are intertwined rather than simply competing, coexisting or favorable [20].

Similar to the cuprate superconductors, in the graphene sheets of the superconducting intercalated graphite CaC${}_6$ [demonstrated in Fig.1(a)], the charge inhomogeneity was also observed [21]. Focusing on every single graphene sheet, the electron density tends to be higher on every third row of hexagons containing a Ca atom [shaded areas in Fig.1(b)], for which we call the stripe region; while the areas separating the stripes are called the domain region. The area of the domain region is twice that of the stripe region. The relationship between superconductivity and the striped charge inhomogeneity in such a system has not yet been investigated theoretically. Quite interestingly, charge inhomogeneity was also observed in the superconducting magic-angle twisted bilayer graphene [22]. So a natural question to ask is, will there be a “universal” relationship between superconductivity and charge inhomogeneity in different superconducting systems?

Fig.1 The schematic plots of (a) the rhombohedral structure of CaC${}_6$ from the scanning tunnelling microscopy (STM) view [21]; (b) the charge distribution on a graphene layer, with higher electron density in the shaded areas; (c) the $d+\mathrm{i}d$-wave pairing factors on the honeycomb lattice; and (d) the full pairing susceptibility $\chi$ and the bare pairing susceptibility $\tilde{\chi }$.

Full size|PPT slide

To answer part of the question, here we address the relationship between superconductivity and charge inhomogeneity in the system of superconducting intercalated graphite CaC${}_6$. Since the striped charge inhomogeneity is hosted by every graphene sheet of CaC${}_6$, we focus on a single graphene sheet, which can be properly described by the Hubbard model on the honeycomb lattice. In this paper we simulate the superconductivity via the quantum Monte Carlo (QMC) calculations of the $d+\mathrm{i}d$-wave pairing, which was suggested to be the dominant pairing symmetry in the Hubbard model on the honeycomb lattice [23, 24]. The charge inhomogeneity is imposed via the external site energy. We consider three types of charge inhomogeneity: Case I corresponds to the experimental finding in Ref. [21], i.e., the total electron density is $\rho =1.2$, and with the introduction of the external site energy, the electron density in the stripe region becomes higher than that in the domain region. Cases II and III are artificial, for the purpose of comparison. For case II, we still set the total electron density $\rho =1.2$, but the electron density in the stripe region is lower than that in the domain region. For case III, the electron density in the stripe region is higher than in the domain region, but we set the total electron density $\rho =0.8$, i.e., switch the experimental electron-rich doping to the hole-rich doping. As we can see later in the results part, case II and case III provide identical pairing information, since the parameter settings of those two cases can be transformed to each other by doing particle-hole transformation, which can also be taken as a consistency check of our QMC calculations. We find there is an optimal charge inhomogeneity for the $d+\mathrm{i}d$-wave pairing in the experiment-related case I; while for cases II and III, the $d+\mathrm{i}d$-wave pairing is monotonically suppressed by the imposed charge inhomogeneity. The spin correlation across the stripe region in the three cases is also examined. Our results provide considerable insight into the relationship of superconductivity and charge inhomogeneity, extend the previous work from cuprate superconductors to the graphitic systems.

2 Model and methodology

We consider the Hubbard model on the honeycomb lattice to describe a single graphene sheet. The charge inhomogeneity is introduced via adding an external site energy $V_0$ on the periodical stripes. The model is given in terms of the Hamiltonian

(1)$\begin{array}{rl}H=& -t\sum _{⟨\mathit{i}\mathbf{,}\mathit{j}⟩\sigma }(c_{\mathit{i}\sigma }^{†}{c}_{\mathit{j}\sigma }^{\phantom{†}}+c_{\mathit{j}\sigma }^{†}{c}_{\mathit{i}\sigma }^{\phantom{†}})+U\sum _{\mathit{i}}{n}_{\mathit{i}\uparrow }{n}_{\mathit{i}\downarrow }\\ & -\mu \sum _{\mathit{i}}(n_{\mathit{i}\uparrow }+n_{\mathit{i}\downarrow })-V_0\sum _{\mathit{i}\in \mathcal{P}}{n}_i,\end{array}$

where $t$ is the fermion hopping amplitude between nearest-neighbor (NN) sites on the honeycomb lattice (here, $\mathit{i}$ and $\mathit{j}$ denote lattice vectors), $U$ denotes an onsite repulsion, and $\mu$ the chemical potential that allows to tune the electron density, denoted $\rho$ in the following. $\mathcal{P}$ denotes the periodical regions where the external site energy $V_0$ is introduced. We work in units of $t=1$ in the following.

In the graphene sheet of the superconductor CaC${}_6$, the stripe region was found on every third row of hexagons containing a Ca atom [21] [Fig.1 (b)]. To allow such stripe pattern to form on the honeycomb lattice with periodical boundary condition, the possible lattice size is $L_1\times L_2$, with $L_1$ parallel to the stripe direction and multiple of $4$, and $L_2$ perpendicular to the stripe direction and multiple of $9$. Fig.1(b) shows a $16\times 9$ lattice. We start from the smallest possible lattice size $4\times 9$, and do the simulations up to lattice with $288$ sites. The electron density in the stripe region was found to be higher than in the domain region [21], so a negative external site energy is required, as indicated by the negative sign of the $V_0$ term in Eq. (1). For the superconductor CaC${}_6$, in the case of complete Ca ionization, $\frac{1}{3}$ electrons will be donated to every C atom [25], but it was estimated around $0.2$ electrons will be donated to every C atom due to the incomplete ionization of the Ca atom [21, 26]. So in our study, we focus on the total electron density $\rho =1.2$. The value of $U$ is chosen to be $3$, which is half of the bandwidth of the Hubbard model on the honeycomb lattice, and is close to the local interaction strength $3.3$ from the first principle calculation of the graphene system [27]. As a comparison, we also consider $U=1.6$ which was estimated to be the effective on-site interactions of graphene (including the effect of the nonlocal Coulomb interactions) [28].

The numerical algorithm used in this paper is the finite-temperature determinant quantum Monte Carlo (DQMC) method [29, 30]. In DQMC, the partition function $Z=\mathrm{T}\mathrm{r}{\mathrm{e}}^{-\beta H}$ can be traced out analytically if the exponential is of quadratic form of fermion creation and annihilation operators, where $\beta$ is the inverse temperature. For that purpose, the quartic interacting term of the Hamiltonian Eq. (1) is transformed into quadratic form using the Hubbard-Stratonovich (HS) transformation [31], with the price of introducing the HS fields. Then the classical Monte Carlo technique is applied to sample the HS fields. For each set of parameters in our simulations, the equal-time measurements are measured at every Monte Carlo sweep, while the imaginary-time quantities are measured every 10th sweep, to reduce the auto-correlation of the data. The chemical potential is tuned to achieve the desired electron density to within an accuracy of $O({10}^{-4})$. The inverse imaginary time $\beta$ is divided into slices $\beta =M\mathrm{\Delta }\tau$ with $\mathrm{\Delta }\tau =0.1$, thus the Trotter error which is proportional to $\mathrm{\Delta }{\tau }^2$ is negligible.

Our purpose is to examine how the stripe pattern affect the dominant $d+\mathrm{i}d$-wave pairing symmetry in the Hubbard model on the honeycomb lattice. Within the DQMC simulations, we can directly access the temperature dependence of the pairing susceptibility for the $d+id$-wave pairing channel,

(2)${\chi }_{d+\mathrm{i}d}={\displaystyle \frac{1}{N_s}\sum _{\mathit{i}\mathbf{,}\mathit{j}}{\int }_0^{\beta }\mathrm{d}\tau ⟨{\mathrm{\Delta }}_{d+\mathrm{i}d,\mathit{i}}^{†}(\tau ){\mathrm{\Delta }}_{d+\mathrm{i}d,\mathit{j}}^{\phantom{†}}(0)⟩,}$

where ${\mathrm{\Delta }}_{d+\mathrm{i}d,\mathit{i}}^{†}(\tau )={\mathrm{e}}^{\tau H}{\mathrm{\Delta }}_{d+\mathrm{i}d,\mathbf{i}}^{†}(0){\mathrm{e}}^{-\tau H}$, $N_s$ is the number of lattice sites, $\beta =1/T$ is the inverse temperature. The corresponding local pairing operators are given as

(3)${\mathrm{\Delta }}_{d+\mathrm{i}d,\mathit{i}}={\displaystyle \frac{1}{\sqrt{N_{d+\mathrm{i}d}}}\sum _l{\mathrm{e}}^{\mathrm{i}(l-1)\frac{2\pi }{3}}(c_{\mathit{i}\uparrow }^{}{c}_{\mathit{i}+{\delta }_l\downarrow }-c_{\mathit{i}\downarrow }^{}{c}_{\mathit{i}+{\delta }_l\uparrow }),}$

where the vectors ${\delta }_l$ ($l=1,2,3$) denote the NN lattice directions from a given lattice site, as shown in Fig.1 (c), $N_{d+\mathrm{i}d}=3$ is the normalization factor.

In order to extract the effective $d+id$-wave pairing susceptibility, we also compute the bare pairing contributions ${\tilde{\chi }}_{d+\mathrm{i}d}$, for which the two-particle terms $⟨c_{\mathit{i}\downarrow }^{†}(\tau )c_{\mathit{j}\downarrow }^{\phantom{†}}(0)c_{\mathit{k}\uparrow }^{†}(\tau )c_{\mathit{l}\uparrow }^{}(0)⟩$ that appear in evaluating the ${\chi }_{d+\mathrm{i}d}$ in Eq. (2) are replaced by the decoupled contributions $⟨c_{\mathit{i}\downarrow }^{†}(\tau )c_{\mathit{j}\downarrow }^{\phantom{†}}(0)⟩⟨c_{\mathit{k}\uparrow }^{†}(\tau )c_{\mathit{l}\uparrow }^{}(0)⟩$ [Fig.1(d)]. The effective pairing susceptibilities are then given as ${\chi }_{d+\mathrm{i}d}^{\text{eff}}=$ ${\chi }_{d+\mathrm{i}d}-{\tilde{\chi }}_{d+\mathrm{i}d}$, and where a positive (negative) value of ${\chi }_{d+\mathrm{i}d}^{\text{eff}}$ signals an enhanced (suppressed) tendency towards $d+\mathrm{i}d$-wave pairing channel [32].

The equal-time spin correlation function $C(\mathit{r})=$ $⟨S_{\mathit{i}+\mathit{r}}^{-}{S}_{\mathit{i}}^{+}⟩$ with $S_{\mathit{i}}^{+}=c_{\mathit{i}\uparrow }^{†}{c}_{\mathit{i}\downarrow }^{\phantom{†}}$ is also measured, to investigate that how the spin correlation is affected by the charge inhomogeneity.

3 Results and discussion

We start from the case corresponding to the experimental finding, i.e., the total electron density is $\rho =1.2$, and the electron density in the stripe region is higher. With the increase of the stripe strength $V_0$, we observe that the electrons are “attracted” to the stripe region, as shown in Fig.2(a) for several lattice sizes and temperatures. For the spatial electron distributions, the finite temperature effect can be ignored for the considered $\beta \ge 8$, while the finite size effect can also be ignored except for the lattice sizes with $L_1=4$. Fig.2(b) presents the corresponding average sign $⟨S⟩$, when $⟨S⟩$ gets close to $0$, the accuracy of the imaginary-time measurements such as the pairing susceptibility cannot be guaranteed, which is the main limitation of the DQMC algorithm [33, 34]. In this paper we constrain our simulations in the region where $⟨S⟩\gtrsim 0.3$, to have a well-controlled statistical accuracy.

Fig.2 (a) The electron densities in the stripe region (${\rho }_{\text{str}}$, solid symbols) and in the domain region (${\rho }_{\text{dom}}$, open symbols), as functions of the stripe strength $V_0$, for different temperatures and lattice sizes. (b) The corresponding average sign $⟨S⟩$ of the DQMC algorithm. The parameters $\rho =1.2$ and ${\rho }_{\text{str}}>{\rho }_{\text{dom}}$ correspond to the experimental finding, and an intermediate interaction strength $U=3$ is chosen. Error bars are smaller than the size of the symbols while not shown.

Full size|PPT slide

In Fig.3, we examine the effect of the charge inhomogeneity on the effective $d+\mathrm{i}d$-wave pairing susceptibility, by tuning the stripe strength $V_0$. Fig.3(a) shows such an effect on a $4\times 9$ lattice for different temperatures. It is clear that, with the increase of $V_0$, the $d+\mathrm{i}d$-wave pairing is firstly enhanced and then gets suppressed. This phenomenon becomes more obvious when the temperature gets lower. The data on a $4\times 36$ lattice are also presented, which show a tiny finite size effect compared with the $4\times 9$ lattice. Then in Fig.3(b), we show the results at $\beta =8$ for more lattice sizes, which exhibit an optimal $V_0\approx 2.5$ for the $d+\mathrm{i}d$-wave pairing. The data on different lattice sizes suggest that the finite size effect of ${\chi }_{d+\mathrm{i}d}^{\text{eff}}$ strongly depends on $L_1$ but is less dependent on $L_2$. The finite size effect is obvious when $L_1$ is small, but for $L_1\ge 16$, the finite size effect is almost negligible, suggesting that the $d+\mathrm{i}d$-wave pairing order will persist in the thermodynamic limit. In total, the effective $d+\mathrm{i}d$-wave pairing susceptibility presented in Fig.3 suggests that the charge distribution pattern observed in the superconductor CaC${}_6$ shows an optimal inhomogeneity for the $d+\mathrm{i}d$-wave pairing, which is the main finding of this work.

Fig.3 The effective $d+\mathrm{i}d$-wave pairing susceptibility as a function of the stripe strength $V_0$, for (a) different temperatures on $4\times 9$ and $4\times 36$ lattices; and (b) $\beta =8$ on different lattice sizes. The parameters $\rho =1.2$ and ${\rho }_{\text{str}}>{\rho }_{\text{dom}}$ correspond to the experimental finding, and an intermediate interaction strength $U=3$ is chosen.

Full size|PPT slide

Such a theoretical finding has the potential to be verified by experiments. For the intercalated graphite CaC${}_6$, one can dope more Ca atoms in the stripe region, and reduce the number of Ca atoms in the domain region. Alternately, in the stripe region, one may replace some of the Ca atoms by the lanthanum (La) atoms; while in the domain region, some of the Ca atoms can be replaced by the potassium (K) atoms. In both ways, more (less) electrons will be donated to the stripe (domain) region. Then one can measure experimentally, to see whether such imposed inhomogeneity can enhance the superconductivity, and whether an optimal inhomogeneity for the superconductivity exists.

To make the conclusion more convincing, we further study two artificial cases: case II with electron-rich doping ($\rho =1.2$), and the stripe strength $V_0$ in Eq. (1) is negative (${\rho }_{\text{str}}<{\rho }_{\text{dom}}$); case III with hole-rich doping ($\rho =0.8$), and the stripe strength $V_0$ is positive (${\rho }_{\text{str}}>{\rho }_{\text{dom}}$). For case II, as shown in Fig.4 for $\beta =8$ and different lattice sizes, the onset of the stripe monotonically suppresses the $d+\mathrm{i}d$-wave pairing, in contrast with case I. The parameter setting of case III can be actually converted into case II by doing the particle-hole symmetric transformation, as indicated by the electron distribution shown in the inset of Fig.4. Indeed, the pairing information we obtain from case III is identical to the case II, as can be seen from the DQMC data of the effective pairing susceptibility (golden empty circles shown in Fig.4). The obvious discrepancy between case I and case II (and case III) strongly suggests that the charge stripe pattern observed in the superconductor CaC${}_6$ [Fig.1(b)] exhibits an optimal inhomogeneity for the $d+\mathrm{i}d$-wave pairing, and such an effect cannot be repeated from artificial charge inhomogeneities.

Fig.4 Effect of stripe strength $|V_0|$ on the effective $d+\mathrm{i}d$-wave pairing susceptibility, for two artificial parameter settings: case II with $\rho =1.2$ and ${\rho }_{\text{str}}<{\rho }_{\text{dom}}$; and case III with $\rho =0.8$ and ${\rho }_{\text{str}}>{\rho }_{\text{dom}}$. The inverse temperature $\beta =8$ and interaction strength $U=3$ are chosen. The data on $4\times 9$ lattice is divided by a factor of $2$, to increase the readability of this plot. The inset shows the evolution of the electron densities in the stripe regions and in the domain regions, for both case II and case III.

Full size|PPT slide

The different effect of case I and case II (and case III) is closely linked to the relationship of the charge/spin density wave order (which is mainly the AF order in the Hubbard model studied in this paper) and superconductivity. Though the relationship of the density wave order and superconductivity was argued to be multifacetted [20], recent studies suggest that the density wave order can enhance superconductivity by serving as an effective pairing “glue” [18, 35, 36]. Considering mainly the domain region (which is twice as large as the stripe region), the increase of $V_0$ in case I drives the electron density in the domain region approach and then get below half-filling [Fig.2(a)]. So there exhibits an optimal $V_0$ for the AF order and also for the superconductivity. While for case II (and case III), the electron density in the domain region is always driven away from half-filling, so the AF order is further destroyed with the increase of $V_0$, and the superconductivity keeps to be suppressed. Such an explanation is certainly rough, but can be taken as a starting point to understand such complicated phenomena.

An interesting feature of the striped phases in the cuprate superconductor described by the Hubbard model on the square lattice is the “$\pi$ phase shift” across the stripe region [37], i.e., the antiferromagnetic (AF) order changes sign as the stripe is traversed. Here we investigate whether such “$\pi$ phase shift” occurs in the Hubbard model on a honeycomb lattice, for all the three cases considered above. According to the stripe pattern observed in the superconductor CaC${}_6$, we measure the equal-time spin correlation function $C(r=4)$ for the lattice sites across the stripe region, as indicated by the green bond in the inset of Fig.5. If the “$\pi$ phase shift” occurs with the onset of the stripe strength $V_0$, $C(r=4)$ would vary from positive value to negative value. However, for the parameters considered here, the signal of the spin correlation becomes pretty weak for $r\ge 3$. Instead, we measure $C(r=1)$ and $C(r=2)$, i.e., the NN and the next-nearest-neighbor (NNN) spin correlation function, as indicated by the blue and pink bonds in the inset of Fig.5(a) and Fig.5, separately. For $r=1$, as we can see in Fig.5(a), the NN spin correlation for case I is suppressed; while for cases II and III, the NN spin correlation is firstly enhanced and then suppressed. For all the three cases, $C(r=1)$ does not change sign. For $r=2$, as we can see in Fig.5(b), the NNN spin correlation for case I is suppressed and changes sign when $V_0$ is strong enough; while for cases II and III, the NNN spin correlation is firstly enhanced and then suppressed, but it doesn’t change sign. As an intuitive explanation, the electron density in the stripe region for case I will be further apart from half-filling [${\rho }_{\text{str}}>1.2$, Fig.2(a)] as $V_0$ increases, hence the AF spin correlation is further destroyed; while for cases II and III, ${\rho }_{\text{str}}$ will firstly approach half-filling and then go away from half-filling (inset of Fig.4), hence the AF spin correlation is firstly enhanced and then suppressed. For case I, when ${\rho }_{\text{str}}$ goes further apart from half-filling, $C(r=1)$ can still keep it’s NN AF order, while $C(r=2)$ begins to lose the longer-range AF order. In this sense, the result suggests that the “$\pi$ phase shift” may occur for the experiment-related case I; while for the artificial cases II and III, no hint of “$\pi$ phase shift” appears.

Fig.5 The spin correlation function $C(r)$ as a function of the stripe strength $|V_0|$ for the three cases considered: case I with $\rho =1.2$ and ${\rho }_{\text{str}}>{\rho }_{\text{dom}}$; case II with $\rho =1.2$ and ${\rho }_{\text{str}}<{\rho }_{\text{dom}}$; and case III with $\rho =0.8$ and ${\rho }_{\text{str}}>{\rho }_{\text{dom}}$, for (a) $r=1$ and (b) $r=2$.

Full size|PPT slide

In the end, we further do the simulation for a different $U=1.6$, which was estimated to be the effective on-site interactions of graphene after including the effect of the nonlocal Coulomb interactions [28]. For this weaker interaction strength $U=1.6$, we can reach lower temperatures. The effective $d+\mathrm{i}d$-wave pairing susceptibilities as functions of the stripe strength are shown in Fig.6. It is clear that for case I, there is an optimal charge inhomogeneity for the superconducting $d+\mathrm{i}d$-wave pairing, which becomes more robust with the lowering of the temperature. In contrast, the charge inhomogeneities of cases II and III monotonically suppress and finally destroy the $d+\mathrm{i}d$-wave pairing.

Fig.6 Effect of stripe strength on the effective $d+\mathrm{i}d$-wave pairing susceptibility, for (a) case I, $U=1.6$, different temperatures on a $4\times 9$ lattice; and (b) all three cases, $U=1.6$, $\beta =12$ on a $4\times 9$ lattice.

Full size|PPT slide

The data shown in Fig.6 for $U=1.6$ are in qualitative accordance with the data shown in Fig.3 for $U=3$. However, it is worth to notice that the optimal values of $V_0$ for the $d+\mathrm{i}d$-wave pairing are different for those two cases, the former ($V_0\approx 1.5$) is smaller than the latter ($V_0\approx 2.5$). This is due to the fact that the $d+\mathrm{i}d$-wave pairing is driven by the interaction strength. For the weak interaction strength ($U=1.6$), one needs some external force ($V_0$) to enhance and reach the optimal $d+\mathrm{i}d$-wave pairing; while for the intermediate interaction strength ($U=3$), one needs a stronger external force to enhance the $d+\mathrm{i}d$-wave pairing equivalently.

4 Conclusions

For the charge inhomogeneity (stripe order) observed in the intercalated graphitic superconductor CaC${}_6$, we find that, with the increase of the stripe strength, the superconducting $d+\mathrm{i}d$-wave pairing is firstly enhanced and then suppressed, and the “$\pi$ phase shift” may occur. As comparisons, we also study two artificial charge inhomogeneities, which tend to suppress the $d+\mathrm{i}d$-wave pairing monotonically, and no signal of the “$\pi$ phase shift” appears.

Our work suggests that there’s an intrinsic universal relationship between the charge inhomogeneity and the superconductivity. Indeed, the charge inhomogeneity, which occurs when several physical interactions — spin, charge, lattice, and/or orbita — are simultaneously active, was observed in many different superconducting systems including cuprate, intercalated graphite and magic-angle twisted bilayer graphene superconductors. Similar to the studies of the cuprate systems based on the Hubbard model on the square lattice [16, 36], we find an “optimal” inhomogeneity for the superconducting pairing in the Hubbard model on the honeycomb lattice, which suggests that the spontaneously formed inhomogeneities in the superconductors play an essential positive role in the mechanism of superconductivity.

References

[1] E. Dagotto . Complexity in strongly correlated electronic systems. Science, 2005, 309(5732): 257

[2] J. Tranquada , B. Sternlieb , J. Axe , Y. Nakamura , S. Uchida . Evidence for stripe correlations of spins and holes in copper oxide superconductors. Nature, 1995, 375(6532): 561

[3] J. M. Tranquada , J. D. Axe , N. Ichikawa , A. R. Moodenbaugh , Y. Nakamura , S. Uchida . Coexistence of, and competition between, superconductivity and charge-stripe order in La_1.6−xNd_0.4Sr_xCuO₄. Phys. Rev. Lett., 1997, 78(2): 338

[4] T. Wu , H. Mayaffre , S. Krämer , M. Horvatić , C. Berthier , W. Hardy , R. Liang , D. Bonn , M. H. Julien . Magnetic-field-induced charge-stripe order in the high-temperature superconductor YBa₂Cu₃O_y. Nature, 2011, 477(7363): 191

[5] G. Ghiringhelli , M. Le Tacon , M. Minola , S. Blanco Canosa , C. Mazzoli , N. Brookes , G. De Luca , A. Frano , D. Hawthorn , F. He , T. Loew , M. M. Sala , D. C. Peets , M. Salluzzo , E. Schierle , R. Sutarto , G. A. Sawatzky , E. Weschke , B. Keimer , L. Braicovich . Long-range incommensurate charge fluctuations in (Y, Nd) Ba₂Cu₃O_6+x. Science, 2012, 337(6096): 821

[6] A. J. Achkar , R. Sutarto , X. Mao , F. He , A. Frano , S. Blanco-Canosa , M. Le Tacon , G. Ghiringhelli , L. Braicovich , M. Minola , M. Moretti Sala , C. Mazzoli , R. Liang , D. A. Bonn , W. N. Hardy , B. Keimer , G. A. Sawatzky , D. G. Hawthorn . Distinct charge orders in the planes and chains of ortho-III-ordered YBa₂Cu₃O_6+δ superconductors identified by resonant elastic X-ray scattering. Phys. Rev. Lett., 2012, 109(16): 167001

[7] J. Chang , E. Blackburn , A. Holmes , N. B. Christensen , J. Larsen , J. Mesot , R. Liang , D. Bonn , W. Hardy , A. Watenphul , M. Zimmermann , E. M. Forgan , S. M. Hayden . Direct observation of competition between superconductivity and charge density wave order in YBa₂Cu₃O_6.67. Nat. Phys., 2012, 8(12): 871

[8] Q. Li , M. Hucker , G. D. Gu , A. M. Tsvelik , J. M. Tranquada . Two-dimensional superconducting fluctuations in stripe-ordered La_1.875Ba_0.125CuO₄. Phys. Rev. Lett., 2007, 99(6): 067001

[9] E. Berg , E. Fradkin , E. A. Kim , S. A. Kivelson , V. Oganesyan , J. M. Tranquada , S. C. Zhang . Dynamical layer decoupling in a stripe-ordered high-T_c superconductor. Phys. Rev. Lett., 2007, 99(12): 127003

[10] T. Valla , A. Fedorov , J. Lee , J. Davis , G. Gu . The ground state of the pseudogap in cuprate superconductors. Science, 2006, 314(5807): 1914

[11] J. M. Tranquada , G. D. Gu , M. Hücker , Q. Jie , H. J. Kang , R. Klingeler , Q. Li , N. Tristan , J. S. Wen , G. Y. Xu , Z. J. Xu , J. Zhou , M. V. Zimmermann . Evidence for unusual superconducting correlations coexisting with stripe order in La_1.875Ba_0.125CuO₄. Phys. Rev. B, 2008, 78(17): 174529

[12] E.ArrigoniS.A. Kivelson, Optimal inhomogeneity for superconductivity, Phys. Rev. B 68, 180503(R) (2003)

[13] I.MartinD.PodolskyS.A. Kivelson, Enhancement of superconductivity by local inhomogeneities, Phys. Rev. B 72, 060502(R) (2005)

[14] K. Aryanpour , E. R. Dagotto , M. Mayr , T. Paiva , W. E. Pickett , R. T. Scalettar . Effect of inhomogeneity on s-wave superconductivity in the attractive Hubbard model. Phys. Rev. B, 2006, 73(10): 104518

[15] E. Berg , E. Fradkin , S. A. Kivelson . Theory of the striped superconductor. Phys. Rev. B, 2009, 79(6): 064515

[16] T. A. Maier , G. Alvarez , M. Summers , T. C. Schulthess . Dynamic cluster quantum Monte Carlo simulations of a two-dimensional Hubbard model with stripelike charge-density-wave modulations: Interplay between inhomogeneities and the superconducting state. Phys. Rev. Lett., 2010, 104(24): 247001

[17] R. Mondaini , T. Ying , T. Paiva , R. T. Scalettar . Determinant quantum Monte Carlo study of the enhancement of d-wave pairing by charge inhomogeneity. Phys. Rev. B, 2012, 86(18): 184506

[18] B. X. Zheng , C. M. Chung , P. Corboz , G. Ehlers , M. P. Qin , R. M. Noack , H. Shi , S. R. White , S. Zhang , G. K. L. Chan . Stripe order in the underdoped region of the two-dimensional Hubbard model. Science, 2017, 358(6367): 1155

[19] B. Ponsioen , S. S. Chung , P. Corboz . Period 4 stripe in the extended two-dimensional Hubbard model. Phys. Rev. B, 2019, 100(19): 195141

[20] E. Fradkin , S. A. Kivelson . High-temperature superconductivity: Ineluctable complexity. Nat. Phys., 2012, 8(12): 864

[21] K. C. Rahnejat , C. A. Howard , N. E. Shuttleworth , S. R. Schofield , K. Iwaya , C. F. Hirjibehedin , C. Renner , G. Aeppli , M. Ellerby . Charge density waves in the graphene sheets of the superconductor CaC₆. Nat. Commun., 2011, 2(1): 558

[22] Y. Jiang , X. Lai , K. Watanabe , T. Taniguchi , K. Haule , J. Mao , E. Andrei . Charge order and broken rotational symmetry in magic-angle twisted bilayer graphene. Nature, 2019, 573(7772): 91

[23] T.MaZ.HuangF.HuH.Q. Lin, Pairing in graphene: A quantum Monte Carlo study, Phys. Rev. B 84, 121410(R) (2011)

[24] T. Ying , S. Wessel . Pairing and chiral spin density wave instabilities on the honeycomb lattice: A comparative quantum Monte Carlo study. Phys. Rev. B, 2018, 97(7): 075127

[25] M. Calandra , F. Mauri . Theoretical explanation of superconductivity in C₆Ca. Phys. Rev. Lett., 2005, 95(23): 237002

[26] T. Valla , J. Camacho , Z. H. Pan , A. V. Fedorov , A. C. Walters , C. A. Howard , M. Ellerby . Anisotropic electron−phonon coupling and dynamical nesting on the graphene sheets in superconducting CaC₆ using angleresolved photoemission spectroscopy. Phys. Rev. Lett., 2009, 102(10): 107007

[27] T. O. Wehling , E. Sasioglu , C. Friedrich , A. I. Lichtenstein , M. I. Katsnelson , S. Blugel . Strength of effective Coulomb interactions in graphene and graphite. Phys. Rev. Lett., 2011, 106(23): 236805

[28] M. Schüler , M. Rosner , T. O. Wehling , A. I. Lichtenstein , M. I. Katsnelson . Optimal Hubbard models for materials with nonlocal Coulomb interactions: Graphene, silicene, and benzene. Phys. Rev. Lett., 2013, 111(3): 036601

[29] R. Blankenbecler , D. J. Scalapino , R. L. Sugar . Monte Carlo calculations of coupled Boson-fermion systems (I). Phys. Rev. D, 1981, 24(8): 2278

[30] S. R. White , D. J. Scalapino , R. L. Sugar , E. Y. Loh , J. E. Gubernatis , R. T. Scalettar . Numerical study of the two-dimensional Hubbard model. Phys. Rev. B, 1989, 40(1): 506

[31] J. E. Hirsch . Two-dimensional Hubbard model: Numerical simulation study. Phys. Rev. B, 1985, 31(7): 4403

[32] S. R. White , D. J. Scalapino , R. L. Sugar , N. E. Bickers , R. T. Scalettar . Attractive and repulsive pairing interaction vertices for the two-dimensional Hubbard model. Phys. Rev. B, 1989, 39(1): 839

[33] E. Y. Loh , J. E. Gubernatis , R. T. Scalettar , S. R. White , D. J. Scalapino , R. L. Sugar . Sign problem in the numerical simulation of many-electron systems. Phys. Rev. B, 1990, 41(13): 9301

[34] M. Troyer , U. J. Wiese . Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations. Phys. Rev. Lett., 2005, 94(17): 170201

[35] E. W. Huang , C. B. Mendl , S. Liu , S. Johnston , H. C. Jiang , B. Moritz , T. P. Devereaux . Numerical evidence of fluctuating stripes in the normal state of high-T_c cuprate superconductors. Science, 2017, 358(6367): 1161

[36] H. C. Jiang , S. A. Kivelson . Stripe order enhanced superconductivity in the Hubbard model. Proc. Natl. Acad. Sci. USA, 2022, 119(1): e2109406119

[37] M. Qin , T. Schafer , S. Andergassen , P. Corboz , E. Gull . The Hubbard model: A computational perspective. Annu. Rev. Condens. Matter Phys., 2022, 13(1): 275

Acknowledgements

We would like to thank R. Mondaini for stimulating discussions. This work was sponsored by the Natural Science Foundation of Chongqing, China (Grant No. cstc2021jcyj-msxmX1009) and the Joint Guiding Project of Natural Science Foundation of Heilongjiang Province (Grant No. LH2019A011).

RIGHTS & PERMISSIONS

2023 Higher Education Press