Unconventional chiral d-wave superconducting state in strained graphene^*

Graphene is a two-dimensional (2D) electronic system on a honeycomb lattice whose electronic excitations behave as massless Dirac quasi-particles. Graphene has been one of the most exciting novel materials in the new century because of its interesting physical properties of massless Dirac fermion and non-trivial structural physics of 2D carbon planes.^[1] A large number of exotic states have been proposed theoretically by including the electron–electron interaction and correlation effects in graphene.^[2] There are also a large number of theoretical works on unconventional superconductivity in graphene, such as p + ip superconducting phase,^[3,4] chiral d-wave superconducting state.^[5–13] Recently, one promising route to induce intrinsic superconductivity in graphene is to reconstruct pseudo Landau levels with the application of strain fields.^[14,15] Flat band superconductivity has been proposed in strained graphene, in which a flat band emerges in the normal state.^[16,17] Flat bands can be seen as a route to high-temperature superconductivity, the critical temperature T_c depends linearly on the electron–phonon coupling constant.^[18–22] Especially, a superconducting state is discovered in a flat band arising in slightly twisted bilayer graphene with strong correlation effects in the recent experiments.^[23,24]

It is widely accepted that both the mechanism and the pairing symmetry in the conventional superconductors are different from high-temperature superconductivity discovered in cuprates.^[25] In high-temperature superconductors, many competing orders, such as charge density wave (CDW), pair density wave (PDW) and spin-density wave, have been proposed to understand the pseudogap phase.^[26] Pair density wave is the Cooper pairs that have a nonzero center-of-mass momentum, and it is always accompanied by CDW.^[27,28] A charge-density wave is a density modulation of electrons, both PDW and CDW can break the translational symmetry at low temperatures.^[29] Studying superconductivity in strained graphene is helpful for understanding unconventional superconductivity in high-temperature superconductors.

In this paper, we explore chiral d-wave superconductivity with CDW and PDW coexisting in doped graphene under strain, in which flat bands emerge in the normal state. It is found that chiral d-wave superconductivity can be stabilized under strain even for slightly doped graphene and its critical temperature is linearly proportional to the strength of the electron–electron interaction. We derive a phase diagram that shows different coexisting states of chiral d-wave superconductivity, CDW and PDW with two different periods in the mean-filed level. It is obvious that these orders may coexist with others and CDW is a subsidiary order. Particularly, we show a superconducting phase only in the presence of CDW and PDW. Some similar conclusions have been given in the previous work,^[17] such as a single PDW state and the coexisting states. However, we emphasize on the superconducting critical temperature and the thermal effect on orders in this paper.

Graphene is made out of carbon atoms arranged in a hexagonal structure, the structure can be seen as a triangular lattice with a basis of two atoms per unit cell. The three nearest-neighbor vectors in real space are give by

$$\begin{eqnarray}&&{\delta }_{1}=\frac{a}{2}(1,\sqrt{3}),{\delta }_{2}=\frac{a}{2}(1,-\sqrt{3}),{\delta }_{3}=-a(1,0),\end{eqnarray} \tag{ 1 }$$

where a = 1.42 Å is the carbon–carbon distance.

We show the structure of graphene in Fig. 1. Flat-band superconductivity is discovered in slightly twisted bilayer graphene, indicating that graphene may host unconventional superconductivity under appropriate conditions.^[23,24] Unlike the twisted bilayer graphene that requires fine tuning of the twisted angle, flat bands can be formed topologically by strain and robustly induced as Landau levels because of the corresponding pseudomagnetic field generated by the strain.^[19] Flat bands in strained graphene have been observed with the strain imposed or engineered by external stretching or periodic ripples in the experiments.^[20,21] We consider the graphene under periodic strain, the strain can be induced by ripple with fixed period L or by external stretching. The tight-binging Hamiltonian for electrons in strained graphene considering that electrons can hop to nearest-neighbor atoms is as follows:

$$\begin{eqnarray}\begin{array}{lll}{H}_{0} & = & -\displaystyle \sum _{i,{j}_{1},\alpha }t{a}_{i,\alpha }^{\dagger }{b}_{{j}_{1},\alpha }-\displaystyle \sum _{i,{j}_{2},\alpha }t{a}_{i,\alpha }^{\dagger }{b}_{{j}_{2},\alpha }\\ & & -\displaystyle \sum _{i,{j}_{3},\alpha }(t+\frac{\beta }{L}\cos \frac{2\pi {x}_{i}}{L}){a}_{i,\alpha }^{\dagger }{b}_{{j}_{3},\alpha }+{\rm{H}}.{\rm{c}}.,\end{array}\end{eqnarray} \tag{ 2 }$$

where a_i,α (${a}_{i,\alpha }^{\dagger }$) annihilates (creates) an electron with spin α on site R_i on sublattice A (an equivalent definition is used for sublattice B, but j = 1,2,3 for three nearest-neighbor sublattices B), t = 2.8 eV is the nearest-neighbor hopping energy, x represents the position coordinates of site R_i, β is a parameter describing the strength of the strain and L = 3n a/2 is the strain period and is conformed to the period of graphene. We use n = 15 in our calculation compared to L = 0.1–10 nm in the experiment and δ t/t ranges from 0 to 0.5.^[20,21,30] It is obvious that there is a new period in the x direction in real space and forms a superlattice, which leads to flat bands in the k_x direction in momentum space. To achieve the flat band state for the strain field experienced by the Dirac electrons, the effective vector potential is A = (0, A_y(x), 0), where ${A}_{y}(x)=\frac{\beta }{L}\cos \frac{2\pi x}{L}$ is the same as the previous work.^[16] We show energy bands of strained graphene in Fig. 2. It is obvious that there are zero energy points in the lowest band in the k_y direction and almost full zero energy bands in the k_x direction. For higher energy bands, there are pseudo Landau levels in the k_y direction and they almost keep constant in the k_x direction. Flat bands are proposed as a route to increase the critical temperature of superconductivity, the superconducting critical temperature T_c is linear with the electron–phonon coupling constant, which is different from the exponential relation for conventional superconductors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To include the correlation effects, the Hamiltonian describing the electron–electron interaction is given by

$$\begin{eqnarray}&&H={H}_{0}+U\displaystyle \sum _{\sigma =i,j}{\hat{n}}_{\sigma \uparrow }{\hat{n}}_{\sigma \downarrow }.\end{eqnarray} \tag{ 3 }$$

In the strong interacting limit when U is large, the Hilbert space is energetically decomposed into singly occupied and doubly occupied spaces so that the electron operator can be decomposed as ${C}_{i,\alpha }^{\dagger }=(1-{n}_{i\bar{\alpha }}){C}_{i,\alpha }^{\dagger }+{n}_{i\bar{\alpha }}{C}_{i,\alpha }^{\dagger }$, where $\bar{\alpha }=-\alpha$. We decompose H as follows:

$$\begin{eqnarray}&&{H}_{0}={\hat{T}}_{0}+{\hat{T}}_{+1}+{\hat{T}}_{-1},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{0} & = & -t\displaystyle \sum _{i,j,\alpha }[(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }^{\dagger }(1-{n}_{j\bar{\alpha }}){b}_{j,\alpha }+{n}_{i\bar{\alpha }}{a}_{i,\alpha }^{\dagger }{n}_{j\bar{\alpha }}{b}_{j,\alpha }\\ & & +(1-{n}_{j\bar{\alpha }}){b}_{j,\alpha }^{\dagger }(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }+{n}_{j\bar{\alpha }}{b}_{j,\alpha }^{\dagger }{n}_{i\bar{\alpha }}{a}_{i,\alpha }],\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{+1} & = & -t\displaystyle \sum _{i,{j}_{1},\alpha }[{n}_{i\bar{\alpha }}{a}_{i,\alpha }^{\dagger }(1-{n}_{j,\bar{\alpha }}){b}_{{j}_{1},\alpha }\\ & & +{n}_{j\bar{\alpha }}{b}_{j,\alpha }^{\dagger }(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }],\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{-1} & = & -t\displaystyle \sum _{i,j,\alpha }[(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }^{\dagger }{n}_{j\bar{\alpha }}{b}_{j,\alpha }\\ & & +(1-{n}_{j\bar{\alpha }}){b}_{j,\alpha }^{\dagger }{n}_{i\bar{\alpha }}{a}_{i,\alpha }].\end{array}\end{eqnarray} \tag{ 7 }$$

If the canonical transformation is generated by S, the transformed Hamiltonian H' will read

$$\begin{eqnarray}&&{H}^{^{\prime} }={{\rm{e}}}^{{\rm{i}}S}H{{\rm{e}}}^{-{\rm{i}}S}=H+[{\rm{i}}S,H]+[{\rm{i}}S,[{\rm{i}}S,H]]+\cdots \end{eqnarray} \tag{ 8 }$$

By requiring ${\hat{T}}_{+1}+{\hat{T}}_{-1}+[{\rm{i}}S,{H}_{U}]=0$, where ${H}_{U}=U\displaystyle \sum _{\sigma =i,j}{\hat{n}}_{\sigma \uparrow }{\hat{n}}_{\sigma \downarrow }$, we have

$$\begin{eqnarray}&&{\rm{i}}S=\frac{1}{U}({\hat{T}}_{+1}-{\hat{T}}_{-1}).\end{eqnarray} \tag{ 9 }$$

Because of no double occupancy, we have

$$\begin{eqnarray}&&{H}^{^{\prime} }={\hat{T}}_{0}+\frac{1}{U}[{\hat{T}}_{+1},{\hat{T}}_{-1}].\end{eqnarray} \tag{ 10 }$$

In the end, we obtain

$$\begin{eqnarray}\begin{array}{lll}{H}^{^{\prime} } & = & {\hat{T}}_{01}+{\hat{T}}_{02}+{\hat{T}}_{03}+\frac{4{t}^{2}}{U}\displaystyle \sum _{i,{j}_{1}}[{S}_{ai}{S}_{b{j}_{1}}\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{1}}^{b}]\\ & & +\frac{4{t}^{2}}{U}\displaystyle \sum _{{}_{i,{j}_{2}}}[{S}_{ai}{S}_{b{j}_{2}}-\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{2}}^{b}]\\ & & +\frac{4}{U}\displaystyle \sum _{{}_{i,{j}_{3}}}{t}_{i{j}_{3}}{t}_{i{j}_{3}}^{* }[{S}_{i}^{a}{S}_{{j}_{3}}^{b}-\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{3}}^{b}],\end{array}\end{eqnarray} \tag{ 11 }$$

where ${t}_{i{j}_{3}}=(t+\frac{\beta }{L}\cos \frac{2\pi {x}_{i}}{L})$ and ${t}_{i{j}_{1}}={t}_{i{j}_{2}}=t$.

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{01} & = & -t\displaystyle \sum _{{}_{i,{j}_{1},\alpha }}[(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }^{\dagger }(1-{n}_{{j}_{1}\bar{\alpha }}){b}_{{j}_{1},\alpha }\\ & & +{n}_{i\bar{\alpha }}{a}_{i,\alpha }^{\dagger }{n}_{{j}_{1}\bar{\alpha }}{b}_{{j}_{1},\alpha }+(1-{n}_{{j}_{1}\bar{\alpha }}){b}_{{j}_{1},\alpha }^{\dagger }(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }\\ & & +{n}_{{j}_{1}\bar{\alpha }}{b}_{{j}_{1},\alpha }^{\dagger }{n}_{i\bar{\alpha }}{a}_{i,\alpha }],\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{02} & = & -t\displaystyle \sum _{{}_{i,{j}_{2},\alpha }}[(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }^{\dagger }(1-{n}_{{j}_{2}\bar{\alpha }}){b}_{{j}_{2},\alpha }\\ & & +{n}_{i\bar{\alpha }}{a}_{i,\alpha }^{\dagger }{n}_{{j}_{2}\bar{\alpha }}{b}_{{j}_{2},\alpha }+(1-{n}_{{j}_{2}\bar{\alpha }}){b}_{{j}_{2},\alpha }^{\dagger }(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }\\ & & +{n}_{{j}_{2}\bar{\alpha }}{b}_{{j}_{2},\alpha }^{\dagger }{n}_{i\bar{\alpha }}{a}_{i,\alpha }],\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}\begin{array}{lll}{\hat{T}}_{03} & = & -\displaystyle \sum _{{}_{i,{j}_{3},\alpha }}[{t}_{i{j}_{3}}(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }^{\dagger }(1-{n}_{{j}_{3}\bar{\alpha }}){b}_{{j}_{3},\alpha }\\ & & +{t}_{i{j}_{3}}{n}_{i\bar{\alpha }}{a}_{i,\alpha }^{\dagger }{n}_{{j}_{3}\bar{\alpha }}{b}_{{j}_{3},\alpha }+{t}_{i{j}_{3}}^{* }(1-{n}_{{j}_{3}\bar{\alpha }}){b}_{{j}_{3},\alpha }^{\dagger }(1-{n}_{i\bar{\alpha }}){a}_{i,\alpha }\\ & & +{t}_{i{j}_{3}}^{* }{n}_{{j}_{3}\bar{\alpha }}{b}_{{j}_{3},\alpha }^{\dagger }{n}_{i\bar{\alpha }}{a}_{i,\alpha }].\end{array}\end{eqnarray} \tag{ 14 }$$

To satisfy the no-double-occupancy constraint, the Guzwiller approximations are adopted by using the renormalized parameters. The low energy Hamiltonian is similar to the t-J model which can be generally expressed as H_eff = H_t + H_J. H_t and H_J can be given by

$$\begin{eqnarray}\begin{array}{lll}{H}_{t} & = & -{t}^{^{\prime} }\displaystyle \sum _{i,{j}_{1},\alpha }{a}_{i,\alpha }^{\dagger }{b}_{{j}_{1},\alpha }-{t}^{^{\prime} }\displaystyle \sum _{{}_{i,{j}_{2},\alpha }}{a}_{i,\alpha }^{\dagger }{b}_{{j}_{2},\alpha }\\ & & -\displaystyle \sum _{{}_{i,{j}_{3},\alpha }}{{t}^{^{\prime} }}_{i{j}_{3}}{a}_{i,\alpha }^{\dagger }{b}_{{j}_{3},\alpha }+{\rm{H}}.{\rm{c}}.,\end{array}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\begin{array}{lll}{H}_{J} & = & \frac{4{t}^{2}}{U}\displaystyle \sum _{{}_{i,{j}_{1}}}[{S}_{ai}{S}_{b{j}_{1}}-\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{1}}^{b}]+\frac{4{t}^{2}}{U}\displaystyle \sum _{{}_{i,{j}_{2}}}[{S}_{ai}{S}_{b{j}_{2}}-\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{2}}^{b}]\\ & & +\frac{4}{U}\displaystyle \sum _{{}_{i,{j}_{3}}}{t}_{i{j}_{3}}{t}_{i{j}_{3}}^{* }[{S}_{i}^{a}{S}_{{j}_{3}}^{b}-\frac{1}{4}{n}_{i}^{a}{n}_{{j}_{3}}^{b}],\end{array}\end{eqnarray} \tag{ 16 }$$

where t' = g_t t, ${{t}^{^{\prime} }}_{3}={g}_{t}(t+\frac{\beta }{L}\cos \frac{2\pi x}{L})$, and ${g}_{t}=\frac{1-n}{1-2{n}_{i\uparrow }{n}_{i\downarrow }/n}=1-n$ describes the level of doping in each carbon atom.

There are different decouple channels for the t-J model for different research projects, we decouple the four fermion terms into the following identity:

$$\begin{eqnarray}\begin{array}{lll}{S}_{i}{S}_{j} & = & -\frac{3}{8}({C}_{i\uparrow }^{\dagger }{C}_{j\uparrow }+{C}_{i\downarrow }^{\dagger }{C}_{j\downarrow })({C}_{j\uparrow }^{\dagger }{C}_{i\uparrow }+{C}_{j\downarrow }^{\dagger }{C}_{i\downarrow })+\frac{3}{8}{n}_{i}\\ & & -\frac{3}{8}({C}_{i\downarrow }^{\dagger }{C}_{j\uparrow }^{\dagger }-{C}_{i\uparrow }^{\dagger }{C}_{j\downarrow }^{\dagger })({C}_{j\uparrow }{C}_{i\downarrow }\\ & & -{C}_{j\downarrow }{C}_{i\uparrow })-\frac{1}{2}{S}_{i}{S}_{j}.\end{array}\end{eqnarray} \tag{ 17 }$$

Lastly, we obtain the useful Hamiltonian for superconductivity in doped graphene as follows:

$$\begin{eqnarray}&&{H}_{sc}=-\displaystyle \sum _{i,j,\alpha }{t}_{i{j}^{^{\prime} }}{a}_{i,\alpha }^{\dagger }{b}_{j,\alpha }-\frac{3{t}_{ij}^{2}}{2U}{\chi }_{ij}^{\dagger }{\chi }_{ij}-\frac{3{t}_{ij}^{2}}{2U}{\Delta }_{ij}^{\dagger }{\Delta }_{ij}.\end{eqnarray} \tag{ 18 }$$

We shall compute the order parameters in a mean approximation with

$$\begin{eqnarray}&&AB\approx \langle A\rangle B+A\langle B\rangle -\langle A\rangle \langle B\rangle .\end{eqnarray} \tag{ 19 }$$

Because of the existence of periodic strain, the mean-field values of the order parameters are set as

$$\begin{eqnarray}\begin{array}{lll}{\chi }_{1} & = & \displaystyle \sum _{{}_{\sigma }}\langle {a}_{i\alpha }^{\dagger }{b}_{{j}_{1}\alpha }\rangle ={\chi }_{10}+{\chi }_{01}{{\rm{e}}}^{{\rm{i}}2\pi x/L}+{\chi }_{01}^{* }{{\rm{e}}}^{-{\rm{i}}2\pi x/L},\\ {\chi }_{2} & = & {\chi }_{1},\\ {\chi }_{3} & = & \displaystyle \sum _{\sigma }\langle {a}_{i\alpha }^{\dagger }{b}_{{j}_{3}\alpha }\rangle ={\chi }_{10}+{\chi }_{03}{{\rm{e}}}^{{\rm{i}}2\pi x/L}+{\chi }_{03}^{* }{{\rm{e}}}^{-{\rm{i}}2\pi x/L},\\ {\Delta }_{1} & = & \langle {b}_{{j}_{1}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{1}\downarrow }{a}_{i\uparrow }\rangle ={\Delta }_{3}{{\rm{e}}}^{{\rm{i}}2\pi /3},\\ {\Delta }_{2} & = & \langle {b}_{{j}_{2}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{2}\downarrow }{a}_{i\uparrow }\rangle ={\Delta }_{3}{{\rm{e}}}^{{\rm{i}}3\pi /4},\\ {\Delta }_{3} & = & \langle {b}_{{j}_{3}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{3}\downarrow }{a}_{i\uparrow }\rangle \\ & = & {\Delta }_{10}+{\Delta }_{01}{{\rm{e}}}^{{\rm{i}}\pi x/L}+{\Delta }_{01}^{* }{{\rm{e}}}^{-{\rm{i}}\pi x/L}\\ & & +{\Delta }_{11}{{\rm{e}}}^{{\rm{i}}2\pi x/L}+{\Delta }_{11}^{* }{{\rm{e}}}^{-{\rm{i}}2\pi x/L},\end{array}\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}\begin{array}{lll}\langle {n}_{ai}\rangle & = & \eta {{\rm{e}}}^{{\rm{i}}2\pi x/L}+{\eta }^{* }{{\rm{e}}}^{-{\rm{i}}2\pi x/L}+{\eta }^{0},\\ \langle {n}_{bj}\rangle & = & \nu {{\rm{e}}}^{{\rm{i}}2\pi x/L}+{\nu }^{* }{{\rm{e}}}^{-{\rm{i}}2\pi x/L}+{\nu }^{0}.\end{array}\end{eqnarray} \tag{ 21 }$$

The above mean values include the chiral d-wave superconducting order, electron hoping term, CDW and PDW with momenta Q = 2π/L and Q/2. The CDW evaluated by |χ₀₁| is the bond field fluctuation in the real space. The CDW has a period of L in real space and Q in k-space, where Q = 2π/L is the CDW ordering wave vector. Δ₁₀ is the uniform chiral d-wave order, Δ₀₁ and Δ₁₁ are the PDW order with period of L and 2L, respectively, corresponding to the momenta of Q and Q/2 of the Cooper pairs. Here η and ν are the doping fluctuations in different units, and they show on-site CDW in the momentum space. In particular, an important approximation has been used in the above assumptions. The rotational symmetry of superconducting order has been broken due to strain. However, we still consider ${\Delta }_{2}={\Delta }_{3}\exp ({\rm{i}}\frac{4\pi }{3})$ and ${\Delta }_{1}={\Delta }_{3}\exp ({\rm{i}}\frac{2\pi }{3})$ in our numerical calculation for simplify. This approximation has effect on the numerical results of superconducting orders although ${\Delta }_{2}={\Delta }_{1}^{* }$ can be retained even under strain, but the qualitative results are unchanged based on this approximation.^[17] The detailed discussions about those orders will be given in the next section. The mean-field Hamiltonian can be written as

$$\begin{eqnarray}\begin{array}{lll}{H}_{{\rm{mean}}} & = & \displaystyle \sum _{{}_{i,j=1,2,3,\alpha }}{t}_{i{j}^{^{\prime} }}{a}_{i,\alpha }^{\dagger }{b}_{j,\alpha }-\frac{3{t}_{ij}^{2}}{2U}{\chi }_{j}({\chi }_{ij}^{\dagger }+{\chi }_{ij})\\ & & -\frac{3{t}_{ij}^{2}}{2U}{\Delta }_{j}({\Delta }_{ij}^{\dagger }+{\Delta }_{ij})+{H}_{{\rm{const}}},\end{array}\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&{H}_{{\rm{const}}}=\displaystyle \sum _{{}_{i,j=1,2,3}}\frac{3{t}_{ij}^{2}}{2U}{\chi }_{j}^{2}+\frac{3{t}_{ij}^{2}}{2U}{\Delta }_{j}^{2}.\end{eqnarray} \tag{ 23 }$$

We note that H_mean is the same as the mean-field Hamiltonian obtained by using the slave-boson method, so both the mean-field theory based on the Gutzwiller approximation used in this paper and the slave-boson method yields the similar results. We perform a discrete Fourier transformation to the mean-field Hamiltonian with

$$\begin{eqnarray}&&\begin{array}{l}{a}_{i\alpha }(r)=\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}{C}_{{A}_{\alpha },{k}_{x},{k}_{y}}{{\rm{e}}}^{{\rm{i}}k\cdot r},\\ {b}_{{j}_{1}\beta }(r+{\delta }_{1})=\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}{C}_{{B}_{\beta },{k}_{x},{k}_{y}}{{\rm{e}}}^{{\rm{i}}k({\delta }_{1}-{\delta }_{3})}{{\rm{e}}}^{{\rm{i}}k\cdot r},\\ {b}_{{j}_{2}\beta }(r+{\delta }_{2})=\displaystyle \sum _{{k}_{x},{k}_{y}}{C}_{{B}_{\beta },{k}_{x},{k}_{y}}{{\rm{e}}}^{{\rm{i}}k\cdot ({\delta }_{2}-{\delta }_{3})}{{\rm{e}}}^{{\rm{i}}k\cdot r},\\ {b}_{{j}_{3}\beta }(r+{\delta }_{3})=\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}{C}_{{B}_{\beta },{k}_{x},{k}_{y}}{{\rm{e}}}^{{\rm{i}}k\cdot r},\end{array}\end{eqnarray} \tag{ 24 }$$

then we obtain

$$\begin{eqnarray}&&{H}_{{\rm{MF}}}=\displaystyle \sum _{k,\alpha }{H}_{t}(k)+{H}_{\Delta }^{{\rm{MF}}}(k).\end{eqnarray} \tag{ 25 }$$

The ground state of the system can be written as

$$\begin{eqnarray}\begin{array}{lll}|\psi (k)\rangle & = & \displaystyle \sum _{j=0}^{2n-1}{W}_{{A}_{\uparrow },j}{C}_{{A}_{\uparrow },{k}_{x}+\frac{\pi j}{L},{k}_{y}}^{\dagger }|0\rangle +{W}_{{B}_{\uparrow },j}{C}_{{B}_{\uparrow },{k}_{x}+\frac{\pi j}{L},{k}_{y}}^{\dagger }|0\rangle \\ & & +{W}_{{A}_{\downarrow },j}{C}_{{A}_{\downarrow },-{k}_{x}+\frac{\pi j}{L},-{k}_{y}}|0\rangle \\ & & +{W}_{{B}_{\downarrow },j}{C}_{{B}_{\downarrow },-{k}_{x}+\frac{\pi j}{L},-{k}_{y}}|0\rangle,\end{array}\end{eqnarray} \tag{ 26 }$$

then the mean-filed Hamiltonian can be generally expressed as

$$\begin{eqnarray}&&{H}_{{\rm{MF}}}=\displaystyle \sum _{{}_{k,\sigma }}{\psi }^{\dagger }(k){h}_{{\rm{MF}}}(k)\psi (k)\end{eqnarray} \tag{ 27 }$$

with h_MF is an 8n × 8n matrix. The optimal values of effective orders can be obtained by the following self-consistent equations:

$$\begin{eqnarray}\begin{array}{lll}{\chi }_{10} & = & \frac{1}{N}\displaystyle \sum _{{}_{i,{j}_{3}}}\langle {a}_{i\uparrow }^{\dagger }{b}_{{j}_{3}\uparrow }+{a}_{i\downarrow }^{\dagger }{b}_{{j}_{3}\downarrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\langle {C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\uparrow },{k}_{x},{k}_{y}}+{C}_{{A}_{\downarrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\downarrow },{k}_{x},{k}_{y}}\rangle,\\ {\chi }_{01} & = & \frac{1}{N}\displaystyle \sum _{{}_{i,{j}_{1}}}\exp (-{\rm{i}}\frac{2\pi }{L}x)\langle {a}_{i\uparrow }^{\dagger }{b}_{{j}_{1}\uparrow }+{a}_{i\downarrow }^{\dagger }{b}_{{j}_{1}\downarrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\exp [{\rm{i}}({k}_{x}+\frac{2\pi }{L})\frac{3a}{2}+i{k}_{y}\frac{\sqrt{3}a}{2}]\\ & & \times \langle {C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\uparrow },{k}_{x},{k}_{y}+\frac{2\pi }{L}}+{C}_{{A}_{\downarrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\downarrow },{k}_{x},{k}_{y}+\frac{2\pi }{L}}\rangle \\ {\chi }_{03} & = & \frac{1}{N}\displaystyle \sum _{{}_{i,{j}_{3}}}\exp (-{\rm{i}}\frac{2\pi }{L}x)\langle {a}_{i\uparrow }^{\dagger }{b}_{{j}_{3}\uparrow }+{a}_{i\downarrow }^{\dagger }{b}_{{j}_{3}\downarrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\langle {C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\uparrow },{k}_{x},{k}_{y}+\frac{2\pi }{L}}+{C}_{{A}_{\downarrow },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\downarrow },{k}_{x},{k}_{y}+\frac{2\pi }{L}}\rangle,\\ {\Delta }_{10} & = & \frac{1}{N}\displaystyle \sum _{{}_{i,{j}_{3}}}\langle {b}_{{j}_{3}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{3}\downarrow }{a}_{i\uparrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\langle {C}_{{B}_{\uparrow },{k}_{x},{k}_{y}}{C}_{{A}_{\downarrow },-{k}_{x},-{k}_{y}}-{C}_{{B}_{\downarrow },-{k}_{x},-{k}_{y}}{C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}\rangle,\\ {\Delta }_{01} & = & \frac{1}{N}\displaystyle \sum _{{}_{i,{j}_{3}}}\exp (-{\rm{i}}\frac{\pi }{L}x)\langle {b}_{{j}_{3}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{3}\downarrow }{a}_{i\uparrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\langle {C}_{{B}_{\uparrow },{k}_{x},{k}_{y}}{C}_{{A}_{\downarrow },-{k}_{x}+\frac{\pi }{L},-{k}_{y}}-{C}_{{B}_{\downarrow },-{k}_{x}+\frac{\pi }{L},-{k}_{y}}{C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}\rangle,\\ {\Delta }_{11} & = & \frac{1}{N}\displaystyle \sum _{i,{j}_{3}}\exp (-{\rm{i}}\frac{2\pi }{L}x)\langle {b}_{{j}_{3}\uparrow }{a}_{i\downarrow }-{b}_{{j}_{3}\downarrow }{a}_{i\uparrow }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y}}}\langle {C}_{{B}_{\uparrow },{k}_{x},{k}_{y}}{C}_{{A}_{\downarrow },-{k}_{x}+\frac{2\pi }{L},-{k}_{y}}\\ & & -{C}_{{B}_{\downarrow },-{k}_{x}+\frac{2\pi }{L},-{k}_{y}}{C}_{{A}_{\uparrow },{k}_{x},{k}_{y}}\rangle,\end{array}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}\begin{array}{lll} & & \frac{1}{N}\displaystyle \sum _{{}_{i,j,\sigma }}\langle {C}_{ai,\sigma }^{\dagger }{C}_{ai,\sigma }+{C}_{bj,\sigma }^{\dagger }{C}_{bj,\sigma }\rangle \\ & = & \frac{1}{N}\displaystyle \sum _{{}_{{k}_{x},{k}_{y},\sigma }}({C}_{{A}_{\sigma },{k}_{x},{k}_{y}}^{\dagger }{C}_{{A}_{\sigma },{k}_{x},{k}_{y}}+{C}_{{B}_{\sigma },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\sigma },{k}_{x},{k}_{y}})\\ & = & {\eta }^{0}+{\nu }^{0},\\ \eta & = & \frac{1}{N}\displaystyle \sum _{{k}_{x},{k}_{y},\sigma }\langle {C}_{{A}_{\sigma },{k}_{x},{k}_{y}}^{\dagger }{C}_{{A}_{\sigma },{k}_{x}+\frac{2\pi }{L},{k}_{y}}\rangle,\\ \nu & = & \frac{1}{N}\displaystyle \sum _{{k}_{x},{k}_{y},\sigma }\langle {C}_{{B}_{\sigma },{k}_{x},{k}_{y}}^{\dagger }{C}_{{B}_{\sigma },{k}_{x}+\frac{2\pi }{L},{k}_{y}}\rangle .\end{array}\end{eqnarray} \tag{ 29 }$$

The temperature dependence of order $\hat{O}$ is obtained by the following equation:

$$\begin{eqnarray}&&O(T)=\frac{1}{N}\displaystyle \sum _{{}_{k}}\langle \hat{O}\rangle =\frac{1}{N}\displaystyle \sum _{{}_{nk}}\langle {\psi }_{nk}|\hat{O}|{\psi }_{nk}\rangle \frac{1}{1+\exp (\frac{{E}_{nk}}{{k}_{{\rm{B}}}T})},\end{eqnarray} \tag{ 30 }$$

where ψ_nk and E_nk are the eigenfunction and eigenvalue of the mean-filed Hamiltonian, k_B is the Boltzmann constant, and T is temperature.

It is natural to extend the d-wave superconducting state in the high-temperature cuprate superconductors with strong Coulomb repulsion to the graphene, due to the sixfold symmetry of the honeycomb lattice, the d-wave superconducting state in the graphene is a spin singlet d_x² − y² ± id_xy-wave state. Some theoretical works have been carried out in the chiral d-wave superconducting state in the doped graphene,^[5–13] but pure graphene is known to be nonsuperconducting until surprising superconductivity is recently discovered in the twisted bilayer graphene. We consider the chiral d-wave superconducting state in the doped graphene with strong electron–electron interaction in the presence of strain. We show our results about the order parameters as a function of temperature in Fig. 3, the coexistence of CDW and PDW must be considered with strain. The CDW and PDW have the same period as strain and the electron population of each carbon atom also has the same period. With the increase of temperature, thermal effect breaks cooper pair, the superconducting order vanishes with PDW. It is obvious that d-wave superconducting state can stably exist even in slight doping δ = 0.08.

Zoom In Zoom Out Reset image size

Compared to the twisted angle in the bilayer graphene, flat energy bands in the graphene can be easily realized in the application of strain. The flat-band superconductivity is a natural extrapolation of the BCS theory. The superconducting critical temperature T_c in conventional superconductors depends exponentially on the product of the electronic density at the Fermi level and the strength of attractive interaction. Flat energy bands enhance the critical temperature which is linearly proportional to the microscopic coupling constant because pairing electrons can abound around a dispersionless energy band. We show critical temperature T_c achieved by Δ(T_c) = 0 with two different strengths of strain versus the microscopic coupling constant g in Fig. 4. It is obvious that the critical temperature linearly depends on the microscopic coupling constant, and the parameter used here (doping = 0.15) is considered as the optimum dopant concentration in the high-temperature superconductors.

Zoom In Zoom Out Reset image size

The pseudogap phase is unclear up to now, but it is considered as a competing phase.^[26] Besides the chiral d-wave superconductivity, many competing orders such as CDW and PDW have been proposed to explain various experimental observations. This is an important characteristic of strongly correlated electron systems. We show the coexistence states of competing orders in Fig. 5. There are chiral d-wave superconductivity, PDW and CDW in the finite temperature limit. Pair density wave with a double period of CDW is non-trivial. It is proved that the emergent Cooper pair order with momentum Q/2 has to be stabilized as the minimum of the free energy.^[17] With the increase of temperature, PDW with momentum Q/2 vanishes firstly, but chiral d-wave superconductivity still exists with CDW. The PDW with momentum Q always accompanies with chiral d-wave superconductivity because of a new period in the momentum space by the strain. According to our calculations, the CDW is an ever-present order. To verify the lowest free energy state, we show the lowest energy band of quasiparticle excitations with complex orders and compared to only present chiral d-wave superconductivity in Fig. 6. It is clear that the system with PDW and CDW has a significant advantage with lower energy, the complex-order system is more stable than the system with only chiral d-wave superconductivity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is a long-sought-after superconducting state with nonvanishing center-of-mass momentum for Cooper pairs, the emergent PDW state is shown to be superconducting and it can be easily realized in our program. We show a superconducting state with only PDW and CDW in Fig. 7, the chiral d-wave superconductivity vanishes with the PDW with momentum Q. The CDW is an auxiliary order, the PDW with momentum Q/2 takes the place of chiral d-wave to become superconducting order. It is proved that nonvanishing superfluid density exists with only PDW. This novel superconducting state is feasible to experimentally realize the graphene.

Zoom In Zoom Out Reset image size

We show different coexisting states of competing orders with a different doping level in Fig. 8. In the slightly doping limit, the chiral d-wave superconductivity vanishes but a superconducting phase in the presence of only CDW and PDW is found. With the increase of doping level, the chiral d-wave superconductivity becomes the main order in the phase diagram. In the case of intermediate doping level, three competing orders coexist in the low-temperature zone, PDW vanishes with the increase of doping. For the doping deeply enough, PDW vanishes even at zero temperature, only the chiral d-wave superconductivity and CDW survive. In summary, there are three different statuses of competing orders: (I) PDW and CDW for slight doping limit, (II) CDW, PDW and chiral d-wave superconductivity, (III) only CDW and chiral d-wave superconductivity.

Zoom In Zoom Out Reset image size

We have investigated unconventional superconducting states in the doped graphene under periodical strain in the strong electron-electron correlated limit. We firstly obtain the effective t–J model with the help of Gutzwiller approximations which satisfies the no-double-occupancy constraint. The optimal values of effective orders are achieved by solving the self-consistent equations in the momentum space. The flat energy bands are included by strain, which is an effective route to increase superconducting critical temperature T_c, which is linearly proportional to the microscopic coupling constant. There are rich phases in the high-temperature superconductor. We show some competing states with different coexistence states of chiral d-wave superconductivity, CDW and PDW. Especially, we show a long-sought-after superconducting state with nonvanishing center-of-mass momentum for Cooper pairs which are well known as the PDW state.

We thank Professor C. Y. Mou and P. H. Chou for fruitful discussions.