Incommensurate charge ordered states in the $\mathit{t-t^{\prime}-J}$ model

We study the incommensurate charge ordered states in the $\mathit{t-t^{\prime}-J}$ model using the Gutzwiller mean field theory on large systems. In particular, we explore the properties of incommensurate charge modulated states referred to as nodal pair density waves (nPDW) in the literature. nPDW states intertwine site and bond charge order with modulated $d$-wave pair order, and are characterized by a nonzero amplitude of uniform pairing; they also manifest a dominant intra-unit cell $d$-density wave form factor. To compare with a recent scanning tunneling microscopy (STM) study [Hamidian \textit{et al.}, Nat. Phys. 3519 (2015)] of the cuprate superconductor BSCCO-2212, we compute the continuum local density of states (LDOS) at a typical STM tip height using the Wannier function based approach. By Fourier transforming Cu and O sub-lattice LDOS we also obtain bias-dependent intra-unit cell form factors and spatial phase difference. We find that in the nPDW state the behavior of form factors and spatial phase difference as a function of energy agrees remarkably well with the experiment. This is in contrast to commensurate charge modulated states, which we show do not agree with experiment. We propose that the nPDW states are good candidates for the charge density wave phase observed in the superconducting state of underdoped cuprates.


I. INTRODUCTION
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 STM 12 , NMR 1 , x-ray diffraction 4 , and resonant x-ray scattering 28 . Charge order appears to be a generic feature of the phase diagram being observed in materials across the cuprate family such as YBCO 2-4 , BSCCO 5 , LaBSCCO 6 , NaCCOC 12 and HBCO 14 . In contrast to the stripe phase observed, e.g. in LBCO 15 , charge order in other cuprates does not generally accompany any magnetic order. These are short range, uni-directional and incommesurate 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. have shown that at lower energies intraunit cell modulations at O sites are in-phase, giving rise to a dominant s -form factor; however, at higher energies, where the signal reaches a peak intensity a d-form factor dominates. The authors of Ref. 13 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 10 . 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 spinfermion model 16,18,19 . Here, charge order appears as a dFF bond order competing with d-wave superconductivity, and having wave vectors along diagonal 16 or horizontal or vertical axes of the Brillouin zone 19 . The PDW state has also been shown to emerge from this model 20 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 21 .
From a more localized perspective, charge ordered states have been investigated in t − J type models using simple mean field approximation, renormalized meanfield theory (Gutzwiller approximation) 11,23 , variational Monte Carlo 22,25,26 and infinite projected pair entangled states methods 24 . 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 11,23,26 . Of these, anti-phase charge density wave (AP-CDW) and nPDW have dominant d-form factor and exist in the doping range where charge order has been experimentally observed 11 .
The AP-CDW is a charge order with commensurate wave vector e.g. [0, .25] or [. 25,0], that has been studied extensively in Ref. 11 and 23. These states have an accompanying superconducting order parameter that forms domains with opposite signs (hence, the name antiphase (AP) charge density wave). The nPDW is an incommensurate charge order with wave vector [0, Q] or [Q, 0], where Q ≈ 0.3. In addition to the modulating component, the pair field has a uniform d-wave component giving rise to nodal structure in density of states at low energies similar to the experimental observation 12,13 . Thus the nPDW intertwines uniform superconductivity, pair density wave 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 9 .
STM and resonant x-ray scattering experiments directly probe the intra-unit cell lattice sites and provide information about form factors 28,29 . However, in most previous theoretical works, the form factor is obtained by finding the renormalization of nearest neighbor hoppings between unit cells along x-and y-axes 16,17,19 . 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 30 . To avoid a three band calculation, and yet be able to compute 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 31,32 .
In this paper, we first self-consistently solve the t − t − J model under Gutzwiller approximation with suitable initial conditions to obtain nPDW state. 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 then calculate the form factors and show that the results for nPDW state are in very good agreement with the STM experiment 9 . We find at low energies that the sform 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.

II. METHOD
The t − t − J Hamiltonian on a 2-D lattice is given by where c † iσ creates an electron at lattice site i = (i x , i y ) with spin σ = ± and S i is the spin operator at this site. P G = i (1 − n i↑ n i↓ ) is the Gutzwiller projection operator, where n iσ = c † iσ c iσ represents the spin dependent number operator for site i. The hopping matrix element t ij is equal to t and t if i and j are nearest and next nearest neighbor sites, respectively. t = −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 34 where the projection operator is replaced by Gutzwiller counting factors leading to the following renormalized Hamiltonian.
(2) For non-magnetic states the simplified Gutzwiller factors are given as 23 The Gutzwiller factors depend on the local values of the hole density δ i , pair field ∆ v ijσ , and bond field χ v ijσ . As in Ref. 11, the superscript v indicates that these quantities are related to, but not same as the physical order parameters. These mean-fields are given as where,σ = −σ. The d-wave superconducting gap order parameter and the bond order along x (y) direction are given as Variational minimization of the ground state energy E g = Ψ 0 |H|Ψ 0 , with respect to the unprojected wavefunction |Ψ 0 under the constraints of normalization and fixed total occupancy leads to the following mean-field Hamiltonian, where Here, δ ij = 1 if i and j are nearest neighbors and 0 otherwise. λ and µ are Lagarange multipliers and N e is the total electron filling.
Since the charge order observed in most of the cuprates is unidirectional 28 , we will focus only on realization of such states in the extended t − 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.
where N is the lattice dimension in y-direction, R iy is the y-component of the lattice vector corresponding to site i, and c ixk creates an electron with transverse momentum k, at the site i x in the 1D lattice. In this "collapsed 1D" representation, the mean-field Hamiltonian becomes where ixjxσ (k) = iy ixiyjx0 e −ikRi y .
A similar expression holds for D ixjxσ (k) and µ ix = µ (ix,0) . The above Hamiltonian can be diagonalized by a spin-generalized Bogoliubov-de Gennes (BdG) transformation, which leads to the following BdG equations, .
(10) Here, ξ ijσ (k) = ijσ (k) − µ iσ δ ij . For simplicity of notation, we have replaced i x , j x with i and j respectively. The mean fields are given following equations, where f represents the Fermi function. Eq. (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 (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 continuum Green's function using 31 .
where W i (r) is the Wannier function at site i, and r is a three dimensional continuum real space vector. The Wannier function employed in this paper was generated using the Wannier90 package 33 , and is similar in form to that used in Ref. 32. Note that the local Green's function contains nonlocal contributions from all lattice sites. The continuum local density of states is now easily obtained as In most of the previous theoretical works 16,17,19 , intraunit cell form factors were calculated using the Fourier transform of the nearest neighbor bond order χ i,i+x(ŷ) , 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 -, and d-form factors as follows.
where F T refers to the Fourier transform and˜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 7 to extract the form factors and found a significant bias dependence 9 . Using the continuum LDOS information, we can follow a similar approach. First, we obtain LDOS Z-maps, defined below, on a plane located at a typical STM tip height (≈ 5Å) above the BiO plane.
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 9,36 , and subsequently assign it to the sublattice Z-maps Cu Z (r, ω), O Z x (r, ω) and O Z y (r, ω). 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 -, and d-form factors as follows.
Another important quantity of interest is the average spatial phase difference (∆φ) between the positive and negative bias energies for the d-form factor modulations. To compute ∆φ in accordance with the experimental procedure 9 , we filter out the characteristic wave vector corresponding to d-form factor modulation (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(r, ω) and determine its phase φ(r, ω). By taking the average of the spatial phase difference at ±ω, we find ∆φ.
whereÕ g x (q, ω) andÕ g y (q, ω) are the Fourier transforms of the sublattice LDOS maps for oxygen x and oxygen y. Width of the Guassian filter was taken to be Λ = 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 ∆ 0 and wave number Q 0 . With this initial guess, the BdG equations (Eqs. 10 and 11) are solved self-consistently until a converged solution is obtained. Then we find the Green's function matrix using Eq. (12) in a supercell set-up to gain higher energy resolution. Finally, we obtain continuum LDOS at height ≈ 5Å above BiO plane using Eqs. (13)- (16), and compute form factors and spatial phase difference using Eqs. (17) and (18), respectively. Fig. 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 − t − J model for the parameters considered. It is wellknown that the Gutzwiller approximation to the homogeneous t − t − 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 vanishes at a critical doping x c ≈ 0.19 which is slightly above the optimal hole doping 8 .

III. RESULTS
It is clear from Fig. 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 11,23,24,26 . We will return to this aspect in section IV, 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 11 where authors find a large number of nearly degenerate inhomogeneous solutions, although they considered a smaller 16 × 16 lattice and t = 0.
In following subsections we describe each charge ordered state in detail.

A. 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 Refs. 11 and 23 (Note that Yang et. al. 23 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 Fig. 2 for a 56 × 56 system at the hole doping x = 0.125. Fig. 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 Fig. 2  The local density of states (LDOS) on lattice sites over a period of APCDW state is plotted in Fig. 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 (FFLO) phase 35 .
In Fig. 2(d) we show the density wave form factors computed from the bond order on nearest neighbor using Eq. 15. Clearly, APCDW state exhibits dominant dform 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.

B. 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 28 . Thus, in order to address experiments in the context of t − t − 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 = 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 Ref. 11, was obtained using a slightly different initialization procedure. In the following we show the results obtained for a 60x60 system with initial wave number guess Q 0 = 0.154. Fig. 3 shows the characteristic features of the nPDW state in real and Fourier space. Fig. 3(a) shows the variation of hole density and gap order parameter in the central region of 60x60 system at a hole doping of 0.125. The root mean square variation in hole density is found to be ∆δ rms = 0.01 holes. Interestingly, an NMR experiment 1 on the charge ordered phase of YBCO, with hole doping of 0.108, finds a charge density variation ∆δ = 0.03±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 order parameter changes sign. The quasi-incommensurate nature of the nPDW is reflected in Fig. 3(b) which shows that the hole density and the order parameter have many Fourier components. The dominant Fourier components in the order parameter (Q ∆ ) and hole density (Q) are 0.15 and 0.3 respectively, satisfying Q = 2Q ∆ as in the APCDW case.

Characterization of nPDW
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, pair density wave, and uniform d-wave superconductivity. Fig. 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 Andreev bound states 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. Fig. 3(d) shows that the form factors computed from bond order parameter have several Fourier components with Q χ = 0.3 the dominant one. Here too, the d-form factor has the largest weight.

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 32 as an input, we compute the continuum LDOS using Eq. 13 and 14. The resulting LDOS map at energy ω = 0.25t and in a 20x20 unit cell area at a height z ≈ 5Å above BiO plane, which is a typical height for STM tip, is shown in Fig. 4(a). Two types of modulating stripe structures can be observed. In Fig.  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.
More importantly, modulations at the two inequivalent O atoms in an unit cell (O x and O y ) are out of phase, i.e. when 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 -form factor. Indeed, a more quantitative analysis of form factors, discussed in following para-graphs, 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 − 0.27t. Remarkably, similar pattern has been observed in the STM experiments 7,9 . In Fig. 4(c) LDOS map is plotted at the negative bias ω = −0.25t in the same region as in (b). Comparing Fig. 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 Ref. [9], this is a characteristic feature of d-form factor density wave. A more quantitative analysis of the phase differences, calculated using Eg. 18, is given in following paragraphs.  [25,25] of a 60 × 60 system located at a height ≈ 5Å 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 Fig. 4(b), and the cell is shown explicitly in the inset of Fig. 5. Similar to the lattice LDOS, two sets of "coherence peaks" at ≈ ±0.21t and 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 23 .

Bias and doping dependence
Also, a small 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 O x and O y atoms, which clearly demonstrates intra-unit cell C 4 symmetry breaking. The difference between the two is maximum at ω ≈ ±0.21t, the scale corresponding to the hybridized Andreev bound states. 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 (Fig. 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 (i0 + ) in Eq. 12 by an energy dependent artificial broadening i0 + +iΓ(ω) where Γ(ω) = α|ω|, as observed in Ref. 38. Fig. 5(b) shows the resulting continuum LDOS spectrum for α = 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 resemble those taken on O x , O y , and Cu sites very closely 12 . We note that the value of spectral gap in BSCCO-2212 spectrum reported in Ref. 12 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  Using the continuum LDOS map, we can calculate energy dependent form factors as formulated in Eq. 17. To calculate the wave vector corresponding to d-form factor modulation (Q d ), we compute the d-form factor (D Z (q, ω)) 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 Q d = [0. 3,0]. This non-dispersing behavior is very similar to that seen in the experiment 9 .
The energy dependence of the form factors at wave vector Q d = [0.3, 0] is now shown in Fig. 6. Similar to the experiment, we find an s -form factor peak at a lower energy and a d-form factor peak at higher energy. Comparing the energy scales in the lattice LDOS (Fig.  3(c)) and continuum LDOS (Fig. 5), we find that the energy at which d-form factor peaks (Ω d ), corresponds to the Andreev bound state peak. By studying the bias dependence of form factors in systems with varying t , 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 Andreev bound state peak in the lattice LDOS. However, the relative weight of the s -and d-form factor depends on the details of band structure and doping. For example if we choose t = 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.
The doping dependence of the peak value of the d-form factor (D Z max ) and corresponding bias (Ω d ) is shown in Fig. 7. Ω d decreases monotonically with hole doping. On the other hand, D Z max shows a non-monotonic behavior as function of hole doping. First, it increases 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 8 .
The average spatial phase difference (∆φ) between the d-form factor density wave modulations at positive and negative biases, computed using Eq. 18, is shown in Fig.  8(a). We find at x = 0.125 that in the vicinity of Fermi level spatial phase difference is zero and turns to π for ω > 0.12t. This bias dependence is in excellent agree- ment with the STM experiment 9 . Fig. 8(b) shows that the energy (Ω π ) at which π phase shift occur decreases with hole doping. We note that in the supplementary information of Ref. 12, the authors show the bias dependence of ∆φ at 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.
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 selfconsistent mean fields in the nPDW state discussed previously. As a first test, we do the following replacements in Eq. 10: δ i → δ 0 and χ v ijσ → χ v 0 , 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 Eq. 10: ∆ v ijσ → ∆ v 0 , 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 Fig. 9(a)-(c), and (d)-(f) respectively. Comparing Fig. 3(c) with Fig. 9(a) and (d), we find that the two sets of coherence peaks in the nPDW state lattice LDOS are indeed originating from the pair density wave, and that the CDW has an insignificant effect. Fig. 9(b) shows that the d-form factor in the pure pair density wave state has the highest magnitude at the energy corresponding to one of the coherence peaks in the lattice LDOS (Ω d = 0.16t) and its bias dependence and overall scale is very similar to the nPDW state (Fig.  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 Fig. 9(e). The importance of PDW is again manifested in Fig. 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 (Fig. 8). However, when the PDW order is set to zero then we get a very different bias dependence of spatial phase difference as evident from Fig. 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.

IV. DISCUSSION
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 − t − J model. However, the energy difference between the uniform state and charge ordered states is of the order of only O(1 meV) per site 11,23 . 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 25,27 . In fact, the short-ranged nature of these states, observed in STM 29 and resonant elastic x-ray scattering experiments 28 , 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 dif- ferent ladder-type domains, as is normally assumed.
As pointed out in Ref. 11, the evolution of the GW 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 10T or 1meV 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 30T may be reflecting this effect 39,40 .
We find that at a given doping, nPDW states with different ordering wave vectors 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 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 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 ) 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 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 2 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 41 , 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 42 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 0 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.

V. CONCLUSIONS
In summary, we have shown that there exist lowenergy, commensurate and incommensurate charge modulated renormalized mean field solutions of the t − t − 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 is 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 NaC-COC. To calculate these patterns, as well as continuum LDOS spectra within the unit cell, we employed the new Wannier function-based method of Ref. 31, 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 -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 with doping in a manner similar to the pseudogap. Note that with the exception of Ref. 11, previous theories of charge ordered states in t − t − 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 RIXS. 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.