Incommensurate charge ordered states in the t–t'–J model

Recent interest in cuprates has been spurred by the observation of charge order in the underdoped regime of its phase diagram, through a variety of experimental tools like  scanning tunneling microscopy (STM) [1], NMR [2], x-ray diffraction [3], and resonant x-ray scattering [4]. Charge order appears to be a generic feature of the phase diagram observed in materials across the cuprate family such as YBCO [3, 5, 6], BSCCO [7], LaBSCCO [8], NaCCOC [1] and HBCO [9]. In contrast to the stripe phase observed, e.g. in LBCO [10], charge order in other cuprates does not generally accompany any magnetic order. These are short range, uni-directional and incommensurate charge modulations with wavevectors $[0,Q]$ or $[Q,0]$ where Q is around 0.3 in reciprocal lattice units.

Understanding the momentum–energy dependence of local intra-unit cell density form factors as a function of doping can give clues to the origin and interplay of charge order with other symmetry broken phases of cuprates. By directly imaging conductance at Cu and O sublattices, Hamidian et al [11] have shown that at lower energies intra-unit cell modulations at O sites are in-phase, giving rise to a dominant $s^{\prime}$-form factor; however, at higher energies, where the signal reaches a peak intensity, a d-form factor dominates. The authors of [12] associated the characteristic energy scale of the d-form factor modulation to the pseudogap scale. Furthermore, in the higher energy range where the d-form factor dominates, the states below and above the Fermi energy show a robust phase difference of π. In a recent development, the same group has also discovered the modulation of the Cooper-pair density in the charge order phase of BSCCO by using Josephson tunneling microscopy [13]. The energy dependence of the form factors, along with the observations of pair density waves (PDW), put constraints on theories of charge order in cuprates.

Theoretical efforts to understand the physics of the charge order phase have mostly focused on the hot-spot approximation or weak-coupling treatment of the spin-fermion model [14–16]. Here, charge order appears as a d-form factor bond order competing with d-wave superconductivity, with wave vectors along diagonal [14] or horizontal or vertical axes of the Brillouin zone [16]. The PDW state has also been shown to emerge from this model [17] as a possible ground state. However, in a strong coupling Eliashberg type treatment of the model, bond order was found to be suppressed for experimentally relevant parameters [18].

From a more localized perspective, charge ordered states have been investigated in $t\mbox{--}J$ type models using renormalized mean-field theory (Gutzwiller approximation) [19, 20], variational Monte Carlo [21–24] and infinite projected entangled pair states methods [25]. In particular, renormalized mean field theory predicts several unidirectional and bidirectional charge ordered states with ground state energies nearly degenerate with the uniform superconducting state [19, 20, 22]. Of these, anti-phase charge density wave (APCDW) and nodal pair density waves (nPDW) have dominant d-form factor and exist in the doping range where charge order has been experimentally observed [20].

The APCDW is a charge order with commensurate wave vector e.g. $[0,0.25]$ or $[0.25,0]$, that has been studied extensively in [19, 20]. These states have an accompanying superconducting order parameter that forms domains with opposite signs (hence, the name anti-phase (AP) charge density wave). The nPDW is an incommensurate charge order with wave vector $[0,Q]$ or $[Q,0]$, where $Q\approx 0.3$. In addition to the modulating component, the pair field has a uniform d-wave component giving rise to nodal structure in the density of states at low energies similar to the experimental observation [1, 12]. Thus the nPDW intertwines uniform superconductivity, PDW and charge order. In this paper, we will show that this state correctly captures the energy dependence of form factors and spatial phase difference as seen in the STM experiments on the superconducting underdoped cuprates [11].

STM and resonant x-ray scattering experiments directly probe the intra-unit cell lattice sites and provide information about form factors [4, 26]. However, in most of previous theoretical works, the form factor is obtained by finding the renormalization of nearest neighbor hoppings between unit cells along x- and y-axes [14, 16, 27]. A three band model approach introducing planar oxygen states can be utilized to resolve this issue, but the treatment of the no-double occupancy criterion is more complicated [28]. To avoid a three band calculation, and yet be able to compute local density of states (LDOS) maps and form factors as in STM, we have utilized a Wannier function based approach that has been successfully used to explain STM conductance maps around impurities in BSCCO and FeSe [29, 30].

In this paper, we first self-consistently solve the $t\mbox{--}{t}^{\prime }\mbox{--}J$ model within the Gutzwiller approximation with suitable initial conditions to obtain APCDW and nPDW states. Next, we compute the lattice Green's functions and then change the basis to continuum space, with Wannier functions as the matrix elements of transformation, to obtain the continuum Green's function. After obtaining the Cu and O sublattice continuum LDOS at typical STM tip heights, we calculate the form factors and show that the results for nPDW state are in very good agreement with the STM experiment [11]. We find at low energies that the ${s}^{\prime }$-form factor has largest contribution; however, at higher energies the d-form factor dominates. Moreover, we show that the d-form factor modulations are in phase at low energies and become out of phase at the energy scale at which d-form factor becomes dominant. Furthermore, we show that the modulation of d-wave pair order is crucial to obtain the bias dependence of the form factors similar to that observed in the experiment.

The $t\mbox{--}t^{\prime} \mbox{--}J$ Hamiltonian on a 2D lattice is given by

$$\begin{eqnarray}&&H=-\displaystyle \sum _{(i,j),\sigma }{P}_{{\rm{G}}}{t}_{{ij}}({c}_{i\sigma }^{\dagger }{c}_{j\sigma }+{\rm{h}}.{\rm{c}}.){P}_{{\rm{G}}}+J\displaystyle \sum _{\langle i,j\rangle }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},\end{eqnarray} \tag{ 1 }$$

where ${c}_{i\sigma }^{\dagger }$ creates an electron at lattice site $i=({i}_{x},{i}_{y})$ with spin $\sigma =\pm$ and ${{\bf{S}}}_{i}$ is the spin operator at this site. ${P}_{{\rm{G}}}={\prod }_{i}(1-{n}_{i\uparrow }{n}_{i\downarrow })$ is the Gutzwiller projection operator, where ${n}_{i\sigma }={c}_{i\sigma }^{\dagger }{c}_{i\sigma }$ represents the spin dependent number operator for site i. $\langle i,j\rangle$ refers to nearest-neighbor sites whereas $(i,j)$ denotes both nearest- and next-nearest-neighbor sites. The hopping matrix element t_ij is equal to t and $t^{\prime}$ if i and j are nearest and next nearest neighbor sites, respectively. $t^{\prime} =-0.3t$ and $J=0.3t$ are used in this work. When necessary to compare with the experimental bias scale, we take t = 400 meV.

The no double occupancy constraint can be treated in Gutzwiller approximation scheme [31] where the projection operator is replaced by Gutzwiller counting factors leading to the following renormalized Hamiltonian.

$$\begin{eqnarray}H & = & -\displaystyle \sum _{(i,j),\sigma }{g}_{{ij}}^{t}{t}_{{ij}}({c}_{i\sigma }^{\dagger }{c}_{j\sigma }+{\rm{h}}.{\rm{c}}.)\\ & & +\displaystyle \sum _{\langle i,j\rangle }\,J\left[{g}_{{ij}}^{s,z}{S}_{i}^{z}{S}_{j}^{z}+{g}_{{ij}}^{s,{xy}}\left(\displaystyle \frac{{S}_{i}^{+}{S}_{j}^{-}+{S}_{i}^{-}{S}_{j}^{+}}{2}\right)\right].\end{eqnarray} \tag{ 2 }$$

For non-magnetic states the simplified Gutzwiller factors are given as [19]

$$\begin{eqnarray}{g}_{{ij}}^{t} & = & {g}_{i}^{t}{g}_{j}^{t},\\ {g}_{i}^{t} & = & \sqrt{\displaystyle \frac{2{\delta }_{i}}{1+{\delta }_{i}}},\\ {g}_{{ij}}^{s,{xy}} & = & {g}_{i}^{s,{xy}}{g}_{j}^{s,{xy}},\\ {g}_{i}^{s,{xy}} & = & \displaystyle \frac{2}{1+{\delta }_{i}},\\ {g}_{{ij}}^{s,z} & = & {g}_{{ij}}^{s,{xy}}={g}_{{ij}}^{s}.\end{eqnarray} \tag{ 3 }$$

The Gutzwiller factors depend on the local values of the hole density ${\delta }_{i}$, pair field ${{\rm{\Delta }}}_{{ij}\sigma }^{v}$, and bond field ${\chi }_{{ij}\sigma }^{v}$. As in [20], the superscript $v$ indicates that these quantities are related to, but not same as the physical order parameters. These mean-fields are given as

$$\begin{eqnarray}{{\rm{\Delta }}}_{{ij}\sigma }^{v} & = & \sigma \langle {{\rm{\Psi }}}_{0}| {c}_{i\sigma }{c}_{j\bar{\sigma }}| {{\rm{\Psi }}}_{0}\rangle ,\\ {\chi }_{{ij}\sigma }^{v} & = & \langle {{\rm{\Psi }}}_{0}| {c}_{i\sigma }^{\dagger }{c}_{j\sigma }| {{\rm{\Psi }}}_{0}\rangle ,\\ {\delta }_{i} & = & 1-\langle {{\rm{\Psi }}}_{0}| {n}_{i}| {{\rm{\Psi }}}_{0}\rangle ,\end{eqnarray} \tag{ 4 }$$

where, $\bar{\sigma }=-\sigma$ and $| {{\rm{\Psi }}}_{0}\rangle$ is the unprojected wavefunction. The d-wave superconducting gap order parameter and the bond order along x (y) direction are given as

$$\begin{eqnarray}{{\rm{\Delta }}}_{i} & = & \displaystyle \frac{1}{8}\displaystyle \sum _{\sigma }({g}_{i,i+\hat{x}}^{t}{{\rm{\Delta }}}_{i,i+\hat{x},\sigma }^{v}+{g}_{i,i-\hat{x},}^{t}{{\rm{\Delta }}}_{i,i-\hat{x},\sigma }^{v}\\ & & -{g}_{i,i+\hat{y}}^{t}{{\rm{\Delta }}}_{i,i+\hat{y},\sigma }^{v}-{g}_{i,i-\hat{y}}^{t}{{\rm{\Delta }}}_{i,i-\hat{y},\sigma }^{v}),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\chi }_{i,i+\hat{x}(\hat{y})}=\displaystyle \frac{1}{2}\displaystyle \sum _{\sigma }{g}_{i,i+\hat{x}(\hat{y})}^{t}{\chi }_{i,i+\hat{x}(\hat{y}),\sigma }^{v}+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ 6 }$$

Variational minimization of the ground state energy ${E}_{{\rm{g}}}=\langle {{\rm{\Psi }}}_{0}| H| {{\rm{\Psi }}}_{0}\rangle$, with respect to the unprojected wavefunction $| {{\rm{\Psi }}}_{0}\rangle$ under the constraints of normalization and fixed total occupancy leads to the following mean-field Hamiltonian,

$$\begin{eqnarray}{H}_{{\rm{MF}}} & = & \displaystyle \sum _{(i,j),\sigma }{\epsilon }_{{ij}\sigma }{c}_{i\sigma }^{\dagger }{c}_{j\sigma }+{\rm{h}}.{\rm{c}}.\\ & & +\displaystyle \sum _{\langle i,j\rangle ,\sigma }\sigma {D}_{{ij}\sigma }^{* }{c}_{i\sigma }{c}_{j\bar{\sigma }}+{\rm{h}}.{\rm{c}}.\,-\,\displaystyle \sum _{i,\sigma }{\mu }_{i\sigma }{n}_{i\sigma },\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray*}{\epsilon }_{{ij}\sigma } & = & -\displaystyle \frac{3}{4}{{Jg}}_{{ij}}^{s}{\chi }_{{ij}\sigma }^{v* }{\delta }_{{ij},\langle {ij}\rangle }-{g}_{{ij}}^{t}{t}_{{ij}}+\displaystyle \frac{\partial W}{\partial {g}_{{ij}}^{s}}\displaystyle \frac{\partial {g}_{{ij}}^{s}}{\partial {\chi }_{{ij}\sigma }^{v}},\\ {D}_{{ij}}^{* } & = & -\displaystyle \frac{3}{4}{{Jg}}_{{ij}}^{s}{{\rm{\Delta }}}_{{ij}\sigma }^{v* }+\displaystyle \frac{\partial W}{\partial {g}_{{ij}}^{s}}\displaystyle \frac{\partial {g}_{{ij}}^{s}}{\partial {{\rm{\Delta }}}_{{ij}\sigma }^{v}},\\ {\mu }_{i\sigma } & = & \mu -2\displaystyle \sum _{j}\left[\displaystyle \frac{\partial W}{\partial {g}_{{ij}}^{s}}\displaystyle \frac{\partial {g}_{{ij}}^{s}}{\partial {n}_{i\sigma }}+\displaystyle \frac{\partial W}{\partial {g}_{{ij}}^{t}}\displaystyle \frac{\partial {g}_{{ij}}^{t}}{\partial {n}_{i\sigma }}\right],\\ W & = & \langle {{\rm{\Psi }}}_{0}| H| {{\rm{\Psi }}}_{0}\rangle -\lambda (\langle {{\rm{\Psi }}}_{0}| {{\rm{\Psi }}}_{0}\rangle -1)-\mu \left(\displaystyle \sum _{i}{n}_{i}-{N}_{{\rm{e}}}\right),\end{eqnarray*}$$

Here, ${\delta }_{{ij},\langle {ij}\rangle }=1$ if i and j are nearest neighbors and 0 otherwise. λ and μ are Lagrange multipliers and N_e is the total electron filling.

Since the charge order observed in most of the cuprates is unidirectional [4], we will focus only on realization of such states in the extended $t\mbox{--}J$ model. We assume that charge modulation is in the x-direction and exploit translational invariance in y-direction by switching to the (${i}_{x},k$) basis defined by following transformation.

$$\begin{eqnarray}&&{c}_{i\sigma }^{\dagger }=\displaystyle \frac{1}{\sqrt{N}}\displaystyle \sum _{k}{c}_{{i}_{x}\sigma }^{\dagger }(k){{\rm{e}}}^{-{{\rm{i}}{kR}}_{{i}_{y}}},\end{eqnarray} \tag{ 8 }$$

where N is the lattice dimension in y-direction, ${R}_{{i}_{y}}$ is the y-component of the lattice vector corresponding to site i, and ${c}_{{i}_{x}\sigma }^{\dagger }(k)$ creates an electron with transverse momentum k and spin σ, at the site i_x in the 1D lattice. In this 'collapsed 1D' representation, the mean-field Hamiltonian becomes

$$\begin{eqnarray}{H}_{{\rm{MF}}} & = & \displaystyle \sum _{\langle {i}_{x},{j}_{x}\rangle ,k,\sigma }{\epsilon }_{{i}_{x}{j}_{x}\sigma }(k){c}_{{i}_{x}\sigma }^{\dagger }(k){c}_{{j}_{x}\sigma }(k)+{\rm{h}}.{\rm{c}}.\\ & & +\displaystyle \sum _{\langle {i}_{x},{j}_{x}\rangle ,k,\sigma }\sigma {D}_{{i}_{x}{j}_{x}\sigma }^{* }(k){c}_{{i}_{x}\sigma }(k){c}_{{j}_{x}\bar{\sigma }}(-k)+{\rm{h}}.{\rm{c}}.\\ & & -\displaystyle \sum _{{i}_{x},k,\sigma }{\mu }_{{i}_{x}\sigma }{n}_{{i}_{x}\sigma }(k),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray*}&&{\epsilon }_{{i}_{x}{j}_{x}\sigma }(k)=\displaystyle \sum _{{i}_{y}}{\epsilon }_{{i}_{x}{i}_{y}{j}_{x}0\sigma }{{\rm{e}}}^{-{{\rm{i}}{kR}}_{{i}_{y}}}.\end{eqnarray*}$$

A similar expression holds for ${D}_{{i}_{x}{j}_{x}\sigma }(k)$ and ${\mu }_{{i}_{x}}={\mu }_{({i}_{x},0)}$. The above Hamiltonian can be diagonalized by a spin-generalized Bogoliubov–de Gennes (BdG) transformation, which leads to the following BdG equations,

$$\begin{eqnarray}\displaystyle \sum _{j}\left[\begin{array}{cc}{\xi }_{{ij}\uparrow }(k) & -{D}_{{ij}\uparrow }(-k)\\ -{D}_{{ji}\uparrow }^{* }(k) & -{\xi }_{{ij}\downarrow }^{* }(-k)\end{array}\right]\left[\begin{array}{c}{u}_{j\uparrow }^{n}(k)\\ {v}_{j\downarrow }^{n}(k)\end{array}\right]={E}_{n\uparrow }(k)\left[\begin{array}{c}{u}_{i\uparrow }^{n}(k)\\ {v}_{i\downarrow }^{n}(k)\end{array}\right].\end{eqnarray} \tag{ 10 }$$

Here, ${\xi }_{{ij}\sigma }(k)={\epsilon }_{{ij}\sigma }(k)-{\mu }_{i\sigma }{\delta }_{{ij}}$. For simplicity of notation, we have replaced i_x, j_x with i and j respectively. The mean fields are given by the following equations,

$$\begin{eqnarray}{n}_{{i}_{x}\uparrow } & = & \displaystyle \frac{1}{N}\displaystyle \sum _{{nk}}| {u}_{{i}_{x}\uparrow }^{n}(k){| }^{2}f({E}_{n\uparrow }(k)),\\ {n}_{{i}_{x}\downarrow } & = & \displaystyle \frac{1}{N}\displaystyle \sum _{{nk}}| {v}_{{i}_{x}\downarrow }^{n}(k){| }^{2}[1-f({E}_{n\uparrow }(-k))],\\ {{\rm{\Delta }}}_{{i}_{x}{j}_{x}\uparrow }^{v}(k) & = & \displaystyle \sum _{n}{u}_{{i}_{x}\uparrow }^{n}(k){v}_{{j}_{x}\downarrow }^{n* }(k)[1-f({E}_{n\uparrow }(k))],\\ {\chi }_{{i}_{x}{j}_{x}\uparrow }^{v}(k) & = & \displaystyle \sum _{n}{u}_{{i}_{x}\uparrow }^{n* }(k){u}_{{j}_{x}\uparrow }^{n}(k)f({E}_{n\uparrow }(k)),\\ {\chi }_{{i}_{x}{j}_{x}\downarrow }^{v}(k) & = & \displaystyle \sum _{n}{v}_{{i}_{x}\downarrow }^{n}(-k){v}_{{j}_{x}\downarrow }^{n* }(-k)[1-f({E}_{n\uparrow }(-k))],\end{eqnarray} \tag{ 11 }$$

where f represents the Fermi function. Equations (10) and (11) are solved self-consistently until the desired accuracy is obtained. With converged values of the eigenfunctions, the Green's function matrix can be calculated using

$$\begin{eqnarray}{G}_{{ij}\sigma }(\omega ) & = & \displaystyle \frac{1}{N}\displaystyle \sum _{k}{g}_{i}^{t}{g}_{j}^{t}{G}_{{i}_{x}{j}_{x}\sigma }(k,\omega ){{\rm{e}}}^{{\rm{i}}{k}({R}_{{i}_{y}}-{R}_{{j}_{y}}),}\\ {G}_{{i}_{x}{j}_{x}\sigma }(k,\omega ) & = & \displaystyle \sum _{n\gt 0}\left[\displaystyle \frac{{u}_{{i}_{x}\sigma }^{n}(k){u}_{{j}_{x}\sigma }^{n* }(k)}{\omega -{E}_{n\sigma }+i{0}^{+}}+\displaystyle \frac{{v}_{{i}_{x}\sigma }^{n* }(k){v}_{{j}_{x}\sigma }^{n}(k)}{\omega +{E}_{n\bar{\sigma }}(k)+i{0}^{+}}\right].\end{eqnarray} \tag{ 12 }$$

Throughout this paper we have used the artificial broadening ${0}^{+}=0.01t$. To compute the LDOS at the STM tip position, we change the basis and obtain the local continuum Green's function using [29].

$$\begin{eqnarray}&&{G}_{\sigma }({\bf{r}},\omega )=\displaystyle \sum _{{ij}}{G}_{{ij}\sigma }(\omega ){W}_{i}({\bf{r}}){W}_{j}^{* }({\bf{r}}),\end{eqnarray} \tag{ 13 }$$

where ${W}_{i}({\bf{r}})$ is the Wannier function at site i, and ${\bf{r}}$ is a three-dimensional continuum real space vector. The Wannier function employed in this paper was generated using the Wannier90 package [32], and is similar in form to that used in [30]. Note that the local Green's function contains nonlocal contributions from all lattice sites. The continuum LDOS is now easily obtained as

$$\begin{eqnarray}&&{\rho }_{\sigma }({\bf{r}},\omega )=-\displaystyle \frac{1}{\pi }{\rm{Im}}[{G}_{\sigma }({\bf{r}},\omega )].\end{eqnarray} \tag{ 14 }$$

In most of the previous theoretical works [14, 16, 27], intra-unit cell form factors were calculated using the Fourier transform of the nearest neighbor bond order ${\chi }_{i,i+\hat{x}(\hat{y})}$, which can be regarded as the measure of charge density at the oxygen atoms on $x(y)$ bonds at lattice site i. We can express s-, ${s}^{\prime }$-, and d-form factors as follows.

$$\begin{eqnarray}{D}^{\chi }({\bf{q}}) & = & {\rm{FT}}({\tilde{\chi }}_{i,i+\hat{x}}-{\tilde{\chi }}_{i,i+\hat{y}})/2,\\ S{^{\prime} }^{\chi }({\bf{q}}) & = & {\rm{FT}}({\tilde{\chi }}_{i,i+\hat{x}}+{\tilde{\chi }}_{i,i+\hat{y}})/2,\\ {S}^{\chi }({\bf{q}}) & = & {\rm{FT}}(1-{\tilde{\delta }}_{i}),\end{eqnarray} \tag{ 15 }$$

where FT refers to the Fourier transform and $\sim$ denotes that the spatial average of the corresponding quantity has been subtracted to emphasize modulating components. Obviously, this quantity does not have any energy dependence. However, STM experiments utilized phase resolved sublattice LDOS information [33] to extract the form factors and found a significant bias dependence [11]. Using the continuum LDOS information, we can follow a similar approach. To analyze this behavior, we first we obtain LDOS Z-maps, defined below, on a plane located at a typical STM tip height ($\approx 5$ Å) above the BiO plane.

$$\begin{eqnarray}&&{\rho }^{Z}({\bf{r}},\omega \gt 0)=\displaystyle \frac{{\displaystyle \sum }_{\sigma }{\rho }_{\sigma }({\bf{r}},\omega )}{{\displaystyle \sum }_{\sigma }{\rho }_{\sigma }({\bf{r}},-\omega )}.\end{eqnarray} \tag{ 16 }$$

Next, we take non-overlapping square regions around each atom in the Z-map, with the size of the region identical to that used in the experiment [11, 34], and subsequently assign it to the sublattice Z-maps ${{\rm{Cu}}}^{Z}({\bf{r}},\omega )$, ${{\rm{O}}}_{x}^{Z}({\bf{r}},\omega )$ and ${{\rm{O}}}_{y}^{Z}({\bf{r}},\omega )$. We note that form factor results are not very sensitive to the size of the square region, however. Here subscripts x and y designate two nonequivalent oxygen atoms in the unit cell in horizontal and vertical directions, respectively. Taking the proper linear combination of the Fourier transform of the sublattice LDOS yields s-, ${s}^{\prime }$-, and d-form factors as follows,

$$\begin{eqnarray}{D}^{Z}({\bf{q}},\omega ) & = & ({\tilde{{\rm{O}}}}_{x}^{Z}({\bf{q}},\omega )-{\tilde{{\rm{O}}}}_{y}^{Z}({\bf{q}},\omega ))/2,\\ S{^{\prime} }^{Z}({\bf{q}},\omega ) & = & ({\tilde{{\rm{O}}}}_{x}^{Z}({\bf{q}},\omega )+{\tilde{{\rm{O}}}}_{y}^{Z}({\bf{q}},\omega ))/2,\\ {S}^{Z}({\bf{q}},\omega ) & = & {\tilde{{\rm{Cu}}}}^{Z}({\bf{q}},\omega ).\end{eqnarray} \tag{ 17 }$$

Another important quantity of interest is the average spatial phase difference (${\rm{\Delta }}\phi$) between the positive and negative bias energies for the d-form factor modulations. To compute ${\rm{\Delta }}\phi$ in accordance with the experimental procedure [9], we filter out the characteristic wave vector corresponding to d-form factor modulation (${{\bf{Q}}}_{d}$) from the continuum LDOS maps at positive and negative energies using a Gaussian filter. Then we take the inverse Fourier transform to obtain the complex spatial map $D({\bf{r}},\omega )$ and determine its phase $\phi ({\bf{r}},\omega )$. By taking the average of the spatial phase difference at $\pm \omega$, we find ${\rm{\Delta }}\phi$.

$$\begin{eqnarray}{D}^{g}({\bf{q}},\omega ) & = & ({\tilde{{\rm{O}}}}_{x}^{g}({\bf{q}},\omega )-{\tilde{{\rm{O}}}}_{y}^{g}({\bf{q}},\omega ))/2,\\ D({\bf{r}},\omega ) & = & \displaystyle \frac{2}{{(2\pi )}^{2}}\displaystyle \int {\rm{d}}{\bf{q}}{{\rm{e}}}^{{\rm{i}}{\bf{qr}}}{D}^{g}({\bf{q}},\omega ){{\rm{e}}}^{-\displaystyle \frac{{({\bf{q}}-{{\bf{Q}}}_{d})}^{2}}{2{{\rm{\Lambda }}}^{2}},}\\ \phi ({\bf{r}},\omega ) & = & \arctan (\mathrm{Im}[D({\bf{r}},\omega )]/\mathrm{Re}[D({\bf{r}},\omega )]),\\ {\rm{\Delta }}\phi & = & \langle \phi ({\bf{r}},\omega )-\phi ({\bf{r}},-\omega )\rangle ,\end{eqnarray} \tag{ 18 }$$

where ${\tilde{{\rm{O}}}}_{x}^{g}({\bf{q}},\omega )$ and ${\tilde{{\rm{O}}}}_{y}^{g}({\bf{q}},\omega )$ are the Fourier transforms of the sublattice LDOS maps for oxygen x and oxygen y. Width of the Guassian filter was taken to be ${\rm{\Lambda }}=1/2N$.

The full computational procedure is summarized as follows. First, we start with trial values for the hole density on each site, bond field on nearest and next nearest neighbor sites and pair field on nearest neighbor sites in a 2D lattice comprising N × N lattice sites. The hole density and bond field are taken to be uniform, whereas the trial pair field is assumed to have a sinusoidal modulation with a given amplitude ${{\rm{\Delta }}}_{0}$ and wave number Q₀. With this initial guess, the BdG equations (equations (10) and (11)) are solved self-consistently until a converged solution is obtained. Then we find the lattice Green's function matrix using equation (12) in a supercell set-up to gain higher energy resolution. Finally, we obtain continuum LDOS at height $\approx \,5$ Å above BiO plane using equations (13) and (14), and compute form factors and spatial phase difference using equations (17) and (18), respectively.

Figure 1 shows the plot of the energy per site (E/t) for the uniform superconducting and the charge ordered states as a function of hole doping. In the inset we plot the gap order parameter in the uniform superconducting state in doping range 0.01–0.48 where it is realized in $t\mbox{--}{t}^{\prime }\mbox{--}J$ model for the parameters considered. It is well-known that the Gutzwiller approximation to the homogeneous $t\mbox{--}t^{\prime} \mbox{--}J$ phase diagram is similar to cuprates, but with a renormalized doping scale. APCDW and nPDW states were found in the doping range of 0.09–0.17 which is below the optimal hole doping of 0.27. We note that in BSCCO charge order is found empirically to vanish at a critical doping ${x}_{{\rm{c}}}\approx 0.19$ which is slightly above the optimal hole doping [35].

Zoom In Zoom Out Reset image size

It is clear from figure 1 that the uniform superconducting state has lower energy per unit site than the charge ordered states. A similar conclusion has been reached in previous studies using variety of numerical techniques [9, 20, 22, 25]. We will return to this aspect in section 4, where we discuss various scenarios in which the energy of charge ordered states can be lowered relative to the uniform superconducting state. Another striking feature is that the APCDW and, nPDW states with different wave vectors are very close in energy at each and every hole doping. This is consistent with the previous study [20] where some of the authors found a large number of nearly degenerate inhomogeneous solutions, although they considered a smaller 16 × 16 lattice and ${t}^{\prime }=0$.

In following subsections, we describe each charge ordered state in detail.

3.1. APCDW

We first present the results for commensurate APCDW state to set the stage for the discussion of the more complicated incommensurate nPDW phase. The APCDW state has been investigated in [19, 20] (note that Yang et al [19] referred to it as the πDW state). To get an APCDW state with a periodicity of four lattice constants, we work with a system size which is multiple of 8 and, initialize the BdG equations (10) and (11) with a pair field modulating at wave number ${Q}_{0}=1/8$. Results are shown in figure 2 for a 56 × 56 system at the hole doping $x=0.125$.

Zoom In Zoom Out Reset image size

Figure 2(a) shows the variation of hole density (δ) and superconducting gap order parameter (Δ) with lattice sites in the central region of the 56 × 56 system. The wavelength of the gap modulation is 8, which is twice of the wavelength of the charge modulation. Moreover, the gap changes sign after each period of the charge modulation, and the hole density is found to be maximum at the domain wall sites where gap order parameter vanishes. As is evident from figure 2(b), these real space observations are reflected in the Fourier domain where we find dominant charge modulation wave vector (Q) and gap modulation wave vector (${Q}_{{\rm{\Delta }}}$) to be 0.25 and 0.125 respectively. Also, the uniform component of the gap is zero. Thus, the APCDW phase intertwines a unidirectional PDW and charge density wave with wave vectors differing by a factor of two. Our result for ${t}^{\prime }=-0.3$ is similar to that of ${t}^{\prime }=0$ in [20].

The LDOS on lattice sites over a period of APCDW state is plotted in figure 2(c) along with the LDOS in homogeneous state. The APCDW LDOS is finite at all lattice sites, and exhibits two sets of coherence peaks with higher energy peak almost coinciding with homogeneous state coherence peaks. The energy peaks are attributed to overlapping Andreev bound states (ABS) which appear at the domain wall sites where superconducting order parameter changes sign [23]. The ABS form a one-dimensional band, and their overlap between neighboring domain walls broadens and shifts the ABS energy away from the chemical potential. Similar conclusions were drawn in connection to the Fulde–Ferrel–Larkin–Ovchinikov phase [36].

In figure 2(d) we show the density wave form factors computed from the bond order on nearest neighbor using equation (15). Clearly, APCDW state exhibits dominant d-form factor at the wave vector $[0.25,0]$. However, the detailed bias dependence of these states do not match experimental results well. We discuss these comparisons in the supplementary information.

3.2. nPDW

It is well established that the charge order observed in most of the cuprates is incommensurate and, generally, does not accompany a magnetic order [4]. Thus, in order to address experiments in the context of $t\mbox{--}{t}^{\prime }\mbox{--}J$ model, we must look for incommensurate order. However, a truly incommensurate charge order can not be realized in a finite lattice calculation. Instead, we can get a quasi-incommensurate state as a mixture of charge ordered states having different commensurate ordering wave vectors. To achieve this, we initialize the BdG equations (10) and (11) with a pair field modulating at wave number which is not a multiple of $1/N$, i.e. ${Q}_{0}\ne \tfrac{m}{N}$ where m is an integer. Such assignment of wave vector ensures that the initial seed to BdG equations has more than one Fourier component, and thus a self-consistent solution may converge to a quasi-incommensurate state. The nPDW state, first reported in [20], was obtained using a slightly different initialization procedure. In the following we show the results obtained for a 60 × 60 system with initial wave number guess ${Q}_{0}=0.154$.

3.2.1. Characterization of nPDW

Figure 3 shows the characteristic features of the nPDW state in real and Fourier space. Figure 3(a) shows the variation of hole density and gap order parameter in the central region of 60 × 60 system at a hole doping of 0.125. The root mean square variation in hole density is found to be ${\rm{\Delta }}{\delta }_{{\rm{rms}}}=0.01$ holes. Interestingly, an NMR experiment [2] on the charge ordered phase of YBCO, with hole doping of 0.108, found a charge density variation ${\rm{\Delta }}\delta =0.03\pm 0.01$ holes which is of the same order as our findings. Like the case of the APCDW here too, the hole density is maximum at the domain wall site which, in this case, is defined as the lattice site where the gap order parameter changes sign. The quasi-incommensurate nature of the nPDW is reflected in figure 3(b) which shows that the hole density and the gap order parameter have many Fourier components. The dominant Fourier components in the order parameter (${Q}_{{\rm{\Delta }}}$) and hole density (Q) are 0.15 and 0.3 respectively, satisfying $Q=2{Q}_{{\rm{\Delta }}}$ as in the APCDW case.

Zoom In Zoom Out Reset image size

Another characteristic feature of the nPDW is the finite uniform component of the gap order parameter which results into non-vanishing d-wave global pairing as opposed to the case of APCDW. Thus, the nPDW state intertwines quasi-incommensurate charge density wave, PDW, and uniform d-wave superconductivity. Figure 3(c) shows the density of states on a number of lattice sites in nPDW state. Like APCDW, there are two sets of coherence peaks. Higher energy peaks are slightly shifted away from the coherence peak in the uniform superconducting state. Lower energy peaks can be attributed to the hybridization of ABS arising due to domains of sign changing pair field. More importantly, the DOS has a nodal structure around Fermi energy at all lattice sites. This is a direct consequence of the non-vanishing global d-wave pairing. Figure 3(d) shows that the form factors computed from bond order parameter have several Fourier components with ${Q}_{\chi }=0.3$ the dominant one. Here too, the d-form factor has the largest weight.

3.2.2. Continuum LDOS

For a more fruitful comparison with STM experiment, we turn to the continuum LDOS and quantities derived from it. With the first-principles Wannier function for BSCCO-2212 [30] as an input, we compute the continuum LDOS using equations (13) and (14). The resulting LDOS map at energy $\omega =0.25t$ and in a 20 × 20 unit cell area at a height $z\approx 5\,\mathring{\rm A}$ above BiO plane, which is a typical height for STM tip, is shown in figure 4(a). Two types of modulating stripe structures can be observed. In figure 4(b) we plot a zoomed-in view of one of these structures in the area bounded by a square as shown in 4(a). Cu and O atoms located in the CuO plane underneath, are represented as dots and open circles, respectively. The LDOS shows modulations around all atoms which, in the Fourier domain, implies that this particular bias has a mixture of all intra-unit cell form factors.

Zoom In Zoom Out Reset image size

More importantly, modulations at the two inequivalent O atoms in an unit cell (${{\rm{O}}}_{x}$ and ${{\rm{O}}}_{y}$) are out of phase, i.e. when ${{\rm{O}}}_{x}$ has large LDOS then O_y has small LDOS. This leads to the conclusion that the d-form factor has larger weight than ${s}^{\prime }$-form factor. Indeed, a more quantitative analysis of form factors, discussed in following paragraphs, shows that the d-form factor has largest weight at this particular bias. This particular pattern is observed in an energy range of $0.21t\mbox{--}0.27t$. Remarkably, a similar pattern has been observed in the STM experiments [11, 33]. In figure 4(c) LDOS map is plotted at the negative bias $\omega =-0.25t$ in the same region as in (b). Comparing figures 4(b) and (c), it is found that the atoms with larger values of LDOS at positive bias have smaller values at negative bias, which implies a spatial phase change of π between positive and negative biases. As emphasized in [11], this is a characteristic feature of d-form factor density wave. A more quantitative analysis of the phase differences, calculated using equation (18), is given in following paragraphs.

3.2.3. Bias and doping dependence

Figure 5(a) shows the bias dependence of continuum LDOS at Cu, ${{\rm{O}}}_{x}$, and ${{\rm{O}}}_{y}$ atomic positions in unit cell [25, 25] of a 60 × 60 system located at a height $\approx \,5\,\mathring{\rm A}$ above the surface BiO plane. The location of this particular unit cell in reference to others can be found in the lower left corner of figure 4(b), and the cell is shown explicitly in the inset of figure 5(a). Similar to the lattice LDOS, two sets of 'coherence peaks' at $\approx \,\pm 0.21t$ and $\pm 0.37t$ can be observed. These peaks correspond to the modulated Andreev state created by the PDW and that associated with the charge density wave energy scale [19].

Zoom In Zoom Out Reset image size

Also, a small V-shaped gap-like feature exists around around the Fermi level due to the uniform component of the gap order parameter. The most striking feature is the difference between the LDOS at ${{\rm{O}}}_{x}$ and ${{\rm{O}}}_{y}$ atoms, which clearly demonstrates intra-unit cell C₄ symmetry breaking. The difference between the two is maximum at $\omega \approx \pm 0.21t$, the scale corresponding to the hybridized ABS. As we will see, this is the bias at which d-form factor has largest magnitude. Another feature of the LDOS is the strong particle-hole asymmetry. Interestingly, this asymmetry is seen to be much more pronounced in the continuum LDOS than lattice LDOS (figure 3(c)).

We expect that these intra-unit cell contrast of these various effects will be mitigated somewhat when nonzero tip size is accounted for [37]. In addition, we expect the higher-energy features to be broadened significantly by inelastic scattering [38]. To see this effect explicitly, we incorporate linear inelastic scattering by replacing the constant artificial broadening term ($i{0}^{+}$) in equation (12) by an energy dependent artificial broadening $i{0}^{+}+i{\rm{\Gamma }}(\omega )$ where ${\rm{\Gamma }}(\omega )=\alpha | \omega |$, as observed in [38]. Figure 5(b) shows the resulting continuum LDOS spectrum for $\alpha =0.25$ . We find that high energy peaks are indeed broadened and can not be resolved any more. This holds for all higher values of α. The spectrum resembles those taken on O_x, O_y, and Cu sites very closely [1]. We note that the value of spectral gap in the BSCCO-2212 spectrum reported in [1] is in the range 80–90 meV for which the value of α is found to be in the range 0.25–0.33, justifying our choice [38].

Using the continuum LDOS map, we can calculate energy dependent form factors as formulated in equation (17). To calculate the wave vector corresponding to d-form factor modulation (${{\bf{Q}}}_{d}$), we compute the d-form factor (${D}^{Z}({\bf{q}},\omega )$) as a function of energy and obtain the wave vector at which it peaks. We find that above a threshold bias, this wave vector does not show any dispersion and remains constant at ${{\bf{Q}}}_{d}=[0.3,0]$. This non-dispersing behavior is very similar to that seen in the experiment [11].

The energy dependence of the form factors at wave vector ${{\bf{Q}}}_{d}=[0.3,0]$ is now shown in figure 6. Similar to the experiment, we find an $s^{\prime}$-form factor peak at a lower energy and a d-form factor peak at higher energy. Comparing the energy scales in the lattice LDOS (figure 3(c)) and continuum LDOS (figure 5), we find that the energy at which d-form factor peaks (${{\rm{\Omega }}}_{d}$), corresponds to the ABS peak. By studying the bias dependence of form factors in systems with varying $t^{\prime}$, doping level and modulation wave vectors, we find that the d-form factor always displays a peak and the particular bias at which it occurs corresponds to the ABS peak in the lattice LDOS. However, the relative weight of the $s^{\prime}$- and d-form factor depends on the details of band structure and doping. For example if we choose $t^{\prime} =0$ at the hole doping 0.125, then d-form factor is found to have largest magnitude at all energies. Lastly, we note that the magnitude of the s-form factor is comparable to others (although it is never the strongest channel), whereas experiment finds it to be smaller than the others.

Zoom In Zoom Out Reset image size

The doping dependence of the peak value of the d-form factor (${D}_{\max }^{Z}$) and corresponding bias (${{\rm{\Omega }}}_{d}$) is shown in figure 7. ${{\rm{\Omega }}}_{d}$ decreases monotonically with hole doping. On the other hand, ${D}_{\max }^{Z}$ shows a non-monotonic behavior as function of hole doping. First, it increases and achieves a maximum at doping $x=0.13$ and then drops rapidly. This is in agreement with the doping dependence of the STM intensity at the density wave modulation wave vector which can be thought as a measure of d-form factor magnitude [35].

Zoom In Zoom Out Reset image size

The average spatial phase difference (${\rm{\Delta }}\phi$) between the d-form factor density wave modulations at positive and negative biases, computed using equation (18), is shown in figure 8(a). We find at $x=0.125$ that in the vicinity of Fermi level spatial phase difference is zero and turns to π for $\omega \gt 0.12t$. This bias dependence is in excellent agreement with the STM experiment [11]. Figure 8(b) shows that the energy (${{\rm{\Omega }}}_{\pi }$) at which π phase shift occur decreases with hole doping. We note that in the supplementary information of [11], the authors show the bias dependence of ${\rm{\Delta }}\phi$ at a few more doping levels, from which one can infer that the energy corresponding to π phase shift decreases with increasing hole doping level, similar to what we observe for the nPDW state.

Zoom In Zoom Out Reset image size

To get a better understanding of the bias dependence of form factors and spatial phase difference, we attempt to disentangle PDW and CDW orders intertwined in the nPDW state, 'by hand'. We start with the self-consistent mean fields in the nPDW state discussed previously. As a first test, we do the following replacements in equation (10): ${\delta }_{i}\to {\delta }_{0}$ and ${\chi }_{{ij}\sigma }^{v}\to {\chi }_{0}^{v}$, where, subscript 0 indicates that the mean fields correspond to the uniform superconducting state. The pair field remains inhomogeneous and unchanged from the nPDW solution. In the second test, we do the following replacements in equation (10): ${{\rm{\Delta }}}_{{ij}\sigma }^{v}\to {{\rm{\Delta }}}_{0}^{v}$, and leave bond field and hole density inhomogeneous and unchanged from the nPDW solution. The chemical potential is adjusted in both tests to yield the correct average electron filling. Results for the lattice LDOS, form factors and spatial phase difference in first and second tests are shown in figures 9(a)–(c), and (d)–(f) respectively. Comparing figure 3(c) with figures 9(a) and (d), we find that the two sets of coherence peaks in the nPDW state lattice LDOS are indeed originating from the PDW, and that the CDW has an insignificant effect. Figure 9(b) shows that the d-form factor in the pure PDW state has the highest magnitude at the energy corresponding to one of the coherence peaks in the lattice LDOS (${{\rm{\Omega }}}_{d}=0.16t$) and its bias dependence and overall scale is very similar to the nPDW state (figure 6). However, when PDW order is artificially set to zero then d-form factor acquires a bias dependence and scale which is very different from the nPDW state as evident from figure 9(e). The importance of the PDW is again manifested in figure 9(c) which shows that setting the charge density modulations to zero artificially has little effect on the spatial phase difference observed in the nPDW state (figure 8(a)). However, when the PDW order is set to zero then we get a very different bias dependence of spatial phase difference as evident from figure 9(f). Thus we conclude that the most significant features in the bias dependence of lattice LDOS, d-form factor and spatial phase difference in the nPDW state are originating from the pair field modulations.

Zoom In Zoom Out Reset image size

Within the inhomogeneous Gutzwiller approximation, for the parameters employed here, the uniform d-wave superconducting state has a lower energy than the charge ordered states at all doping levels. Thus the nPDW is not the ground state of the $t\mbox{--}t^{\prime} \mbox{--}J$ model. However, the energy difference between the uniform state and charge ordered states is of the order of only ${ \mathcal O }$(1 meV) per site [19, 20]. Thus, it is entirely plausible that other effects not included in the model such as disorder and electron–phonon interactions may stabilize these fluctuating charge ordered states [23, 39]. In fact, the short-ranged nature of these states, observed in STM [26] and resonant elastic x-ray scattering experiments [4], suggests that disorder might be playing an important role. Different local disorder environments may then also pin slightly different states, resulting in slightly different local LDOS patterns that can be identified in STM images, not just two different ladder-type domains, as is normally assumed. As pointed out in [20], the evolution of the Gutzwiller factors with doping is responsible for the remarkable degeneracy of the various charge states shown in figure 1 across the doping range. The energy splitting of these states above the homogeneous superconducting state remains almost the same across this range as well. Thus the addition of a magnetic field on the order of 10 T or 1 meV per site can potentially stabilize long-range charge order. It is tempting to conclude that the recent observation of charge order in YBCO with a large correlation length, at a magnetic field of order 30 T may be reflecting this effect [40, 41].

We find that at a given doping, nPDW states with different ordering wave vectors ${\bf{Q}}$ around $[0.3,0]$ exist. Keeping the same initial guess and changing the system size (N × N) results in charge ordered states with slightly different ${\bf{Q}}=[Q,0],$ since Q is a multiple of $1/N$. However, LDOS, form factor and spatial phase difference results are insensitive to such small changes. All such states at nearby ${\bf{Q}}$ are extremely close in energy, and hence, at the level of the Gutzwiller approximation, we can not quantitatively address the doping dependence of the charge order wave vector. However, the bias dependence of the form factors and spatial phase difference is robust with respect to the change of ordering wave vector, band structure ($t^{\prime}$) and doping. We always find a dominant d-form factor at higher energies and a shift of π in the average spatial phase difference beyond a particular energy scale.

The analysis presented in the previous section, whereby PDW and CDW order were artificially suppressed independently, strongly suggests that PDW character is necessary to explain the spectral characteristics, in particular the bias dependence of the intra unit cell form factors and spatial phase difference in experimental measurements on BSCCO.

It is important to note further that the bias dependence of the form factors in the current theory is the clear result of electronic correlations in the CuO₂ plane. It has been observed in x-ray spectroscopy that the plane in YBCO, for example, buckles in a pattern of O displacements that mimics a d-wave form factor [42], and suggested that this structural pattern imprints itself on the local tunneling conductance. However, it is difficult to see how such a structural effect should be sensitive to the applied bias, as seen in experiment and predicted here. Nor is it clear why, in such a scenario, the other form factors can be stabilized in other bias ranges.

Very recently, a new paper [43] appeared analyzing the conductance maps in BSCCO according to a new type of 'phase resolved electronic structure visualization', concluding that the charge states observed were locally commensurate with lattice constant 4a₀ separated by phase slips (discommensurate), rather than incommensurate. The authors also stated that their findings were consistent with strong-coupling r-space based theories rather than Fermi surface driven instabilities. The latter conclusion is consistent with a t–J model description of the underdoped cuprates, and is quite consistent with ours. Since we have not included disorder or allowed for discommensuration, we cannot address their findings directly in this work. However it seems intuitive that the introduction of disorder may favor discommensuration. We leave this for a later project.

In summary, we have shown that there exist low-energy, commensurate and incommensurate charge modulated renormalized mean field solutions of the $t\mbox{--}t^{\prime} \mbox{--}J$ model that are not the ground state at any filling, but which are extremely close to the energy of the uniform superconducting state. Furthermore, the incommensurate charge ordered states, called nPDW are intertwined with modulated superconductivity, and display properties remarkably similar to STM observations of the 1D modulated states seen on the surface of BSCCO and NaCCOC. These are well-established features of cuprate physics that have intrigued workers in the field for almost a decade, but until now have defied explanation. Among these properties are the same spectra and pattern of tunneling conductance maps within the unit cell as observed by STM on under- to optimally doped BSCCO and NaCCOC. To calculate these patterns, as well as continuum LDOS spectra within the unit cell, we employed the new Wannier function-based method of [29], which enables the calculation of the wavefunctions in the correlated state at any 3D position, including several Å above the surface where the STM tip is placed. This gives us an unprecedented ability to compare with details of the experiments in the charge ordered regime.

In addition, the bias dependence of intra-unit cell d-, $s^{\prime}$- and s-form factors and their spatial phase difference were obtained in the nPDW state and display good agreement with the STM observations. The energy of the peak d-wave form factor depends on doping in a manner similar to the pseudogap. Note that with the exception of [20], previous theories of charge ordered states in $t\mbox{--}t^{\prime} \mbox{--}J$ type models treated only commensurate (4 unit cell wavelength) charge order states, and could express observables only in terms of bias-independent bond variables.

We have furthermore argued that, while the incommensurate states found here are not the ground state of the system studied, relatively small perturbations can stabilize them. In particular, we discussed the possibility that impurities stabilize the charge order, leading to the disordered 1D patterns observed in STM on BSCCO and NaCCOC. This disordered ground state would also be consistent with the short-range charge-order observed by resonant x-ray scattering [4]. In such a system, a magnetic field should suppress superconductivity and eventually favor long-range charge order, as observed in experiments. We leave investigation of this intriguing scenario to a future study.

The authors are grateful for useful discussions with B M Andersen, J C Davis, K Fujita, S Hayden, A Kreisel, S Maiti, C Pépin, T M Rice, S Sachdev, A-M S Tremblay and S Verret. PC and PJH were supported by NSF grant NSF-DMR-1407502. PJH's work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293. WLT and TKL were supported by Taiwan Ministry of Science and Technology Grant 104-2112-M-001-005. Part of calculation was supported by the National Center for High Performance Computing in Taiwan.