\section{Introduction}

The pseudogap (PG) phase in cuprate superconductors provides one of
the most active subjects on high-$T_c$ superconductivity.
The PG phase is characterized by highly anomalous properties\cite{timusk99,norman05} 
which are rather universal for all cuprate superconductors. 
One of the puzzling observations 
comes from angle-resolved photoemission spectroscopy (ARPES) 
measurements,\cite{damascelli03} which revealed arc-shaped disconnected 
Fermi surfaces,\cite{norman98}  called Fermi arcs, instead of a large Fermi surface.
In underdoped cuprates the PG opens below a temperature $T^*$, 
which is far above the superconducting transition temperature $T_{sc}$. 
Furthermore, in contrast to the behavior of $T_{sc}$, 
$T^*$ increases with decreasing doping in the underdoped region. 
The PG is very anisotropic along the Fermi surface. 
It has a maximal gap in the $(0,0)$-$(0,\pi)$ direction (antinodal direction)
and vanishes upon approaching the Brillouin zone diagonal (nodal direction), 
similar to the  $d$-wave superconducting gap. 

In spite of the consensus on the existence of the PG, its origin and
nature remain elusive. 
There are two major scenarios. 
One is  that the PG originates from
preformed pairs above $T_{sc}$.\cite{emery95,norman07}  
The other is that the PG is distinct from the
superconducting gap and associated with a certain order which competes 
with superconductivity, but 
both coexist at low temperature, leading to ``two gaps'' in 
the electronic spectrum.\cite{vishik10,kondo11,yoshida12} 
Several phenomenological models which are in favor of the two-gap scenario 
were already studied in various contexts, but invoking different orders,  
such as $d$-wave charge density wave ($d$CDW),\cite{chakravarty01} 
$d$-wave Fermi surface deformations,\cite{yamase09}  
charge density wave\cite{castellani95,becca96,hashimoto10,he11} 
including stripes,\cite{kivelson03,vojta09} 
phase separation (PS),\cite{emery93,castellani95,becca96} 
and others such as 
resonating-valence-bond-type charge order\cite{yang06}  
and loop-current order.\cite{varma0699} 

The $d$CDW is a flux phase, where 
orbital currents flow around each plaquette in a staggered pattern. 
The electronic spectrum in the flux phase has a gap with $d$-wave 
symmetry, the same as the superconducting gap symmetry. 
The flux phase was obtained in the large-$N$ approach to the $t$-$J$ model 
in various formalisms\cite{affleck88a,affleck89,morse91,cappelluti99,foussats04} 
and the presence of flux correlations was confirmed by the exact diagonalization.\cite{leung00}  
On the other hand, in the Hubbard model, 
the dynamical cluster approximation failed to detect 
static long-range order of the $d$CDW\cite{macridin04}  
whereas the variational cluster approximation 
showed that the $d$CDW is a metastable solution.\cite{lu11}  
Fluctuations associated with the $d$CDW 
can provide a route to address the PG. 
A perturbative analysis of the electron self-energy due to $d$CDW fluctuations 
catches many important features observed 
by ARPES, not only a PG and its associated Fermi arcs\cite{greco09,bejas11} but also 
a semiquantitative aspect of renormalization of 
the electron band dispersion in the PG phase.\cite{greco11} 

The $d$-wave Fermi surface deformations are driven by a 
$d$-wave Pomeranchuk\cite{pomeranchuk59} 
instability ($d$PI), leading to an electronic nematic state.\cite{kivelson98} 
In this state, an orientational symmetry of the systems is broken 
without breaking however translational invariance. The $d$PI was found 
in the slave-boson mean-field,\cite{yamase00} exact diagonalization,\cite{miyanaga06} 
and variational Monte Carlo\cite{edegger06} studies in the $t$-$J$ model, 
and also in the Hubbard model.\cite{metzner00,wegner02,buenemann12} 
The $d$PI itself does not become the leading instability in most of theoretical studies. 
However, it was pointed out that the models retain appreciable correlations of 
the $d$PI,\cite{yamase00,okamoto10} 
which then may lead to a giant response to a small $xy$ anisotropy. 
Such a giant response was actually observed in the PG region in 
YBa$_2$Cu$_3$O$_y$, which has a small anisotropy originating from 
the orthorhombic crystal structure, 
by neutron scattering\cite{hinkov07,hinkov08} and transport measurements.\cite{daou10} 
Theoretical studies for the former\cite{yamase06,yamase09} 
and the latter\cite{hackl09} confirmed that idea. 


Charge-stripe order is extensively discussed 
for cuprates.\cite{kivelson03,vojta09} 
Since the charge-stripe order breaks both orientational and translational 
symmetry of the system, the stripe phase is also called an electronic 
smectic phase\cite{kivelson98} and has lower symmetry than the nematic phase.  
The experimental observation of charge order  
in La-based cuprates\cite{tranquada95} provides grounds to consider 
the stripe order. 
A charge-stripe solution was indeed obtained in the 
density matrix renormalization group (DMRG) study\cite{white98,white00} 
in the $t$-$J$ model. However in the presence of 
the second nearest-neighbor hopping integral the charge-stripe order turned out 
to be unstable in the $t$-$J$ model.\cite{white99,tohyama99} 


PS is also another possible instability in the $t$-$J$ model.\cite{marder90,emery90}  
It is however still highly debated whether the model indeed shows the instability 
toward PS\cite{hellberg00,ivanov04} or not\cite{yokoyama96,white98,white00,hellberg99,shih98,putikka00,lugas06} 
in a parameter region realistic to cuprates. 
Although PS is in general suppressed by long-range Coulomb forces, 
strong charge fluctuations in the proximity to PS can be important 
and responsible for anomalous properties in the PG 
phase and superconductivity.\cite{emery93} 
In fact, the proximity to PS plays an important role 
to generate a singular interaction between electrons at zero momentum transfer 
as shown in the infinite-$U$ Hubbard Holstein model.\cite{castellani95,becca96} 
When long-range Coulomb interactions are added, 
the singularity shifts to 
a finite momentum transfer, leading to an incommensurate charge density wave  
similar to stripes.\cite{castellani95,becca96} 


Theoretically it is believed that the 
two-dimensional (2D) $t$-$J$ and Hubbard models  
contain the main ingredients for describing cuprates,\cite{anderson87}  
i.e., antiferromagnetism at
zero doping, a metallic state at finite doping, and a strong tendency to
$d$-wave superconductivity. 
Given that various charge instabilities are invoked to address 
the PG, and also other anomalous properties in cuprates, 
it is interesting to study 
what kind of charge instabilities are favored in the 2D $t$-$J$ model 
by treating all possibilities on equal footing in a controllable scheme. 

In this paper, we analyze the 2D $t$-$J$ model in terms of Hubbard
operators by including the nearest-neighbor Coulomb interaction $V$ 
to avoid a subtle feature of PS; our main results are not 
affected by the presence of $V$. 
We apply a large-$N$ expansion formulated in a path integral 
representation.\cite{foussats00,foussats02,foussats04}  
In this approach the two spin components are
extended to $N$ and an expansion in powers of
the small parameter $1/N$ is performed, 
providing a controllable scheme 
without a perturbative expansion in any model parameter.  
In addition, different kind of instabilities can be studied on equal footing, 
allowing us to perform a stability analysis on
all possible charge instabilities already at leading order. 
We find that the $t$-$J$ model shows tendencies to the flux 
and electronic nematic state in a wide doping region. 
In particular, the nematic state has a strong tendency to become inhomogeneous. 
Close to half-filling, bond-order phase (BOP) and PS are also obtained. 

In the next section, we first provide a brief summary of our theoretical scheme 
and then explain the most important charge instabilities, 
$d$CDW, $d$PI, BOP, and PS. Our results are presented in Sec.~III 
and are discussed by comparing with literature in Sec.~IV. 
Implications for cuprates are also discussed in the same section. 

\section{Theoretical framework}
\subsection{Large-\boldmath{$N$} approach to the \boldmath{$t$-$J$-$V$} model}
In a previous paper,\cite{foussats04}  
a large-$N$ expansion for the $t$-$J$-$V$ model was formulated in terms of 
a path integral representation for the Hubbard $X$-operators. 
For the sake of a self-contained presentation, 
we first summarize the formalism. 

The $t$-$J$-$V$ model is described by the following Hamiltonian, 
\begin{equation}\label{Hc}
H = -\sum_{i, j,\sigma} t_{i\,j}\tilde{c}^\dag_{i\sigma}
\tilde{c}_{j\sigma}
+ J \sum_{\langle i,\, j\rangle} (\vec{S}_i \cdot \vec{S}_j-\frac{1}{4} n_i n_j)
+ V \sum_{\langle i,\, j\rangle} n_i n_j \, ,
\end{equation}
where $t_{i\,j} = t$ $(t')$ is the hopping integral between the 
first (second) nearest-neighbor sites on a square lattice;  
$J$ and $V$ are the exchange interaction
and the Coulomb repulsion, respectively, between the nearest-neighbor sites. 
The main role of the $V$-term in the present study is to suppress the tendency toward 
PS 
while, in other works,\cite{zeyher98,gazza99} the $V$-term was included to investigate 
its effect on superconductivity. 
$\tilde{c}^\dag_{i\sigma}$ and $\tilde{c}_{i\sigma}$ are 
the creation and annihilation operators of electrons 
with spin $\sigma$ ($\sigma=\downarrow$,$\uparrow$),  respectively, 
under the constraint that 
the double occupancy of electrons is excluded 
at any site $i$. $n_i$ is the electron density operator. 

The electron and spin operators are connected to Hubbard operators\cite{hubbard63}  
via  $\tilde{c}^\dag_{i \sigma}=X_i^{\sigma 0}$, 
$\tilde{c}_{i \sigma}=X_i^{0 \sigma}$, 
$S_i^+=X_i^{\uparrow \downarrow}$, 
$S_i^-=X_i^{\downarrow \uparrow}$, 
and $n_i=X_i^{\uparrow \uparrow}+X_i^{\downarrow \downarrow}$. 
The operators $X_i^{\sigma 0}$ and $X_i^{0 \sigma}$ are called fermionlike, 
whereas the operators $X_i^{\sigma \sigma'}$ and
$X_i^{00}$ are called bosonlike; $X_i^{00}$ will be introduced later [\eq{eq:v1}]. 
After writing Hamiltonian (1) in terms of the Hubbard operators,  
we extend the spin degree of freedom to $N$ channels and obtain 
the Hamiltonian in the large-$N$ formalism, 
\begin{equation}
H = - \frac{1}{N}\sum_{i, j, p}\; t_{i j} X_{i}^{p 0}X_{j}^{0p} 
+ \frac{J}{2N} \sum_{\langle i,j\rangle, pp'} (X_{i}^{p p'}X_{j}^{p' p} - X_{i}^{p p} X_{j}^{p' p'}) 
+ \frac{V}{N}\sum_{\langle i,j\rangle, p p'} X_{i}^{p p} X_{j}^{p' p'}
-\mu\sum_{i,p}\;X_{i}^{p p} \label{eq:H} \,.
\end{equation}
The spin index $\sigma$ is extended to a new index $p$, which 
runs from $1$ to $N$. In order to obtain a finite theory in the $N$-infinite limit,  
$t$, $t'$, $J$ and $V$ are rescaled as
$t/N$, $t'/N$, $J/N$ and $V/N$, respectively. 
The chemical potential $\mu$ is introduced in \eq{eq:H}. 

In the path integral formulation our Euclidean Lagrangian reads 
\begin{eqnarray}
L_E =  \frac{1}{2}
\sum_{i, p}\frac{({\dot{X_{i}}}^{0 p}\;X_{i}^{p 0}
+ {\dot{X_{i}}}^{p 0}\;
X_{i}^{0 p})} {X_{i}^{0 0}} + H
\label{lagrangian}
\end{eqnarray}
with the following two additional constraints, 
\begin{eqnarray}
X_{i}^{0 0} + \sum_{p} X_{i}^{p p} - \frac{N}{2}=0 \; , \label{eq:v1}
\end{eqnarray}
and
\begin{eqnarray}
X_{i}^{p p'} - \frac{X_{i}^{p 0} X_{i}^{0p'}}{X_{i}^{0 0}}=
0 \;, \label{eq:v2}
\end{eqnarray}
which are imposed on the path integral via two $\delta$-functions.
In \eq{lagrangian}, ${\dot{X_{i}}}^{p 0}=\partial_\tau {X_{i}}^{p 0}$ 
and $\tau$ is the euclidean time, namely $\tau={\rm i}t$. 
Equation~(\ref{eq:v1}) is the $N$-extended completeness condition.  
The form of the kinetic term in the Lagrangian $L_E$, as well as 
the constraint \eq{eq:v2}, comes from the requirement that  
the $X$-operators should fulfill their commutation rules. 
For details we refer to 
Ref.~\onlinecite{foussats00}. 
In the path integral approach we associate Grassmann and usual
bosonic variables with fermionlike and bosonlike $X$-operators,
respectively. 

We now discuss the main steps needed to introduce a large-$N$
expansion.\cite{foussats02,foussats04} 
First the $V$-term in the Hamiltonian 
is written in terms of $X_{i}^{0 0}$ by using \eq{eq:v1}. 
We then eliminate the bosonic variables $X^{p p'}$ 
by implementing the $\delta$-function associated with \eq{eq:v2}. 
The completeness condition [\eq{eq:v1}] is imposed by 
introducing Lagrange multipliers $\lambda_i$. 
We write $X_{i}^{0 0}$ and $\lambda_i$  in terms of
static mean-field values, $r_0$ and $\lambda_0$, and fluctuation fields, 
$\delta R_i$ and $\delta \lambda_i$, 
\begin{eqnarray}
X_{i}^{0 0} &=& N r_{0}(1 + \delta R_{i}) \nonumber \\
\lambda_{i} &=&\lambda_{0}+ \delta{\lambda_{i}} \,.
\label{lambdai}
\end{eqnarray}
In addition, we introduce the following fermion fields\cite{misc-fandX} defined by 
\begin{eqnarray}\label{fdag}
f^{\dag}_{i p} &=& \frac{1}{\sqrt{N r_{0}}}X_{i}^{p 0} \nonumber \,,\\
f_{i p} &=& \frac{1}{\sqrt{N r_{0}}}\;X_{i}^{0 p}\,.
\end{eqnarray}
The exchange interaction is then described by four fermion fields, which 
are decoupled through a Hubbard-Stratonovich transformation 
by introducing a field associated with a bond variable, 
\begin{equation}
\Delta_{ij} = J \sum_{p} \frac{f^{\dag}_{j p} f_{i p}}{ \sqrt{(1 + \delta R_{i})
(1 + \delta R_{j})}} \,. 
\end{equation}
The field $\Delta_{ij}$ is parameterized by 
\begin{equation}
\Delta_i^{\eta}=\Delta(1+r_i^\eta+iA_i^\eta)\,,
\label{staticDelta}
\end{equation}  
where $r_i^{\eta}$ and $A_i^{\eta}$ correspond to the real and imaginary parts of the 
fluctuations of the bond variable,
respectively, and $\Delta$ is a static mean-field value. 
The index $\eta$ takes two values associated with the bond directions
 ${\eta}_{1}=(1,0)$ and ${\eta}_{2}=(0,1)$ on a square lattice. 
After expanding $1/(1+\delta R)$ in powers of $\delta R$, 
we obtain an effective Lagrangian, which can be written in terms of 
a six-component boson field 
\begin{equation}
\delta X^{a} = (\delta
R\;,\;\delta{\lambda},\; r^{{\eta}_{1}},\;r^{{\eta}_{2}}
,\; A^{{\eta}_{1}},\;
A^{{\eta}_{2}})\,,
\end{equation} 
the fermions $f_{p}$, and their interactions. 


From the quadratic part for fermions we obtain an electronic 
propagator in the paramagnetic phase, 
\begin{eqnarray}\label{G0}
G({\bf k}, {\rm i}\nu_{n}) = \frac{1}{{\rm i}\nu_{n} -\varepsilon_{{\bf k}}}\,.
\end{eqnarray}
Here ${{\bf k}}$ and ${\rm i}\nu_{n}$ are the momentum and fermionic  Matsubara frequency, respectively, and the electronic dispersion $\varepsilon_{{\bf k}}$ is 
\begin{eqnarray}\label{Ek}
\varepsilon_{{\bf k}} = -2(t r_0+\Delta) (\cos k_{x}+\cos k_{y})-
4t' r_0 \cos k_{x} \cos k_{y} - \mu \,. 
\end{eqnarray}
Here $\lambda_0$ in \eq{lambdai} was absorbed in the chemical potential $\mu$. 


From the completeness condition [\eq{eq:v1}] $r_0$ is equal to $\delta/2$,  
where $\delta$ is the hole doping rate away from half-filling. 
The field $\Delta$ is given  
by the expression 
\begin{eqnarray}{\label {Delta}}
\Delta = \frac{J}{4N_s} \sum_{{\bf k}, \eta} \cos(k_\eta) n_F(\varepsilon_{\bf k}) \; , 
\end{eqnarray}
where $n_F$ is the Fermi function and $N_s$ is the total 
number of lattice sites. For a given doping, 
$\mu$ and $\Delta$ are determined self-consistently by solving \eq{Delta} and 
\begin{equation}
(1-\delta)=\frac{2}{N_s} \sum_{{\bf k}} n_F(\varepsilon_{\bf k})\,.
\end{equation}

The quadratic part for $\delta X^{a}$ 
defines a $6\times6$ bare bosonic propagator 
$D^{(0)}_{ab}({\bf q},\mathrm{i}\omega_n)$, which 
after Fourier transformation reads, 
\begin{widetext}
\begin{equation} \label{D0inverse}
[D^{(0)}_{ab}({\bf q},\mathrm{i}\omega_n)]^{-1} = N 
\left(
\begin{array}{llllll}
\frac{\delta^2}{2}\left(V-\frac{J}{2}\right)
[\cos(q_x) + \cos(q_y)] & \frac{\delta}{2} & 0 & 0 & 0 & 0\\
\frac{\delta}{2} & 0 & 0 & 0 & 0 & 0\\
0 & 0 & \frac{4\Delta^2}{J} & 0 & 0 & 0\\
0 & 0 & 0 & \frac{4\Delta^2}{J} & 0 & 0\\
0 & 0 & 0 & 0 & \frac{4\Delta^2}{J} & 0\\
0 & 0 & 0 & 0 & 0 & \frac{4\Delta^2}{J}
\end{array}
\right) \; ,
\end{equation}
\end{widetext}
where ${{\bf q}}$ and ${\rm i}\omega_{n}$ are the momentum and bosonic 
Matsubara frequency, respectively. 
The quantity $D^{(0)}_{ab}({\bf q},\mathrm{i}\omega_n)$ 
describes all possible types of bare charge susceptibilities. 
The bare susceptibilities are renormalized already at leading order 
to become dressed ones,  which are given by the Dyson equation 
\begin{equation}
D^{-1}_{ab}({\bf q},\mathrm{i}\omega_n)
= [D^{(0)}_{ab}({\bf q},\mathrm{i}\omega_n)]^{-1} - \Pi_{ab}({\bf q},\mathrm{i}\omega_n)\,.
\label{dyson}
\end{equation}
Following the diagrammatic rules in Ref.~\onlinecite{foussats04}, 
the $6\times6$ boson self-energies $\Pi_{ab}$ are computed as 
\begin{widetext}
\begin{eqnarray}\label{eq:Piab}
&& \Pi_{ab}({\bf q},\mathrm{i}\omega_n)
            = -\frac{N}{N_s}\sum_{{\bf k}} h_a({\bf k},{\bf q},\varepsilon_{\bf k}-\varepsilon_{{\bf k}-{\bf q}}) 
            \frac{n_F(\varepsilon_{{\bf k}-{\bf q}})-n_F(\varepsilon_{\bf k})}
                                  {\mathrm{i}\omega_n-\varepsilon_{\bf k}+\varepsilon_{{\bf k}-{\bf q}}} 
            h_b({\bf k},{\bf q},\varepsilon_{\bf k}-\varepsilon_{{\bf k}-{\bf q}}) \nonumber \\
&& \hspace{25mm} - \delta_{a\,1} \delta_{b\,1} \frac{N}{N_s}
                                       \sum_{\bf k} \frac{\varepsilon_{\bf k}-\varepsilon_{{\bf k}-{\bf q}}}{2}n_F(\varepsilon_{\bf k}) \; .
\end{eqnarray}
\end{widetext}
The prefactor $N$ in front of the right hand side of Eq.~(17) comes from 
the sum over the $N$ channels of $p$. 
Thus, the $6\times6$ boson self-energies $\Pi_{ab}$ are 
of the same order as $[D^{(0)}_{ab}]^{-1}$ [see Eq.~(15)]. 
In Eq.~(17) $h_a$ is an effective six-component interaction vertex which comes 
from the interaction terms between bosonic and fermionic fields 
derived from the effective Lagrangian. 
The explicit expression for $h_a$  is given by 
\begin{widetext}
\begin{align}
 h_a({\bf k},{\bf q},\nu) =& \left\{
                   \frac{2\varepsilon_{{\bf k}-{\bf q}}+\nu+2\mu}{2}+
                   2\Delta \left[ \cos\left(k_x-\frac{q_x}{2}\right)\cos\left(\frac{q_x}{2}\right) +
                                  \cos\left(k_y-\frac{q_y}{2}\right)\cos\left(\frac{q_y}{2}\right) \right];1;
                 \right. \nonumber \\
               & \left. -2\Delta \cos\left(k_x-\frac{q_x}{2}\right); -2\Delta \cos\left(k_y-\frac{q_y}{2}\right);
                         2\Delta \sin\left(k_x-\frac{q_x}{2}\right);  2\Delta \sin\left(k_y-\frac{q_y}{2}\right)
                 \right\} \; . 
\end{align}
\end{widetext}


From the $N$-extended 
completeness condition [\eq{eq:v1}] we see that the charge operator 
$X^{00}$ is $O(N)$, while the operators $X^{pp}$ are $O(1)$. 
Consequently, the $1/N$ approach emphasizes   
the effective charge interactions. 
In fact, while in leading order  
charge susceptibilities  contain collective effects, they 
enter the spin susceptibilities in the next-to-leading order. 
Similarly, superconductivity appears 
in the next-to-leading order.\cite{cappelluti99,zeyher98}   
Therefore, instabilities of the paramagnetic phase 
are expected only, in leading order, in the charge sector. 

In leading order, our formalism agrees with the $1/N$ slave-boson 
formalism.\cite{morse91}  
However, in the present approach the fermion variables $f_{ip}$ are proportional to the
$X$-operators [\eq{fdag}] and should not be confused with
the spinons in the slave-boson approach. In addition, $\delta R$ [\eq{lambdai}] 
is proportional to charge fluctuations and not related to holons. 
Since the $X$-operators are treated as fundamental objects, 
problems associated with fluctuations of gauge fields 
in the slave-boson approach\cite{lee06} are avoided. 
Our formalism was also checked to yield results consistent with 
the exact diagonalization\cite{merino03,bejas06} 
as well as results in another formalism of 
the $1/N$ expansion in leading order.\cite{cappelluti99}  


\subsection{Instabilities of the paramagnetic phase} 
An instability of the paramagnetic phase is signaled by 
the divergence of the static susceptibilities defined 
by $D_{ab}({\bf q},\mathrm{i}\omega_n)$ for a continuous phase transition. 
Therefore we study eigenvalues and eigenvectors 
of the matrix $[D_{ab}({\bf q}, \rm{i}\omega_n)]^{-1}$ at $\rm{i}\omega_n=0$. 
When an eigenvalue crosses zero at a given doping rate, temperature $T$, 
and momentum ${\bf q}$, an instability occurs toward a phase characterized 
by the corresponding eigenvector.
We have found five instabilities associated with eigenvectors $V^a$ 
explained below.   

\begin{figure}[htp]
\centering
\includegraphics[width=8cm]{figure1.eps}
\caption{Sketch of various phases appearing in our work. 
Commensurate orders for (a) $d$CDW, (b) $d$PI, 
(c) BOP$_x$, and (d) BOP$_{xy}$ in real space. 
The commensurate $d$PI has a momentum ${\bf q}=(0,0)$ 
whereas the commensurate $d$CDW, BOP$_x$, and BOP$_{xy}$ 
have ${\bf q}=(\pi,\pi)$. Solid and dashed lines in (b)-(d) 
represent the strong and weak bonds, respectively.
(e) Fermi surface deformations (black line) 
associated with the commensurate $d$PI; the original Fermi surface 
is sketched by gray lines. 
}
\label{fig1}
\end{figure}


a) $V^a =\frac{1}{\sqrt{2}} (0,0,0,0,1,-1)$, which corresponds to the freeze of  
the imaginary parts of the bond variable [\eq{staticDelta}]. 
The pure imaginary contribution to the hopping term generates 
a net magnetic flux in each plaquette, leading to 
the instability toward the flux or $d$CDW phase as already 
found previously.\cite{affleck88a,affleck89,morse91,cappelluti99,foussats04}
The commensurate flux phase is characterized by the 
modulation vector ${\bf q}=(\pi,\,\pi)$ 
and describes staggered circulating currents as sketched in \fig{fig1} (a),    
whereas the incommensurate phase is characterized by ${\bf q}\ne(\pi,\,\pi)$. 

b) $V^a = \frac{1}{\sqrt{2}}(0,0,1,-1,0,0)$, which corresponds to the freeze of 
the real  parts of the bond variable. 
This eigenvector corresponds to the commensurate [${\bf q}=(0\,,0)$] or 
incommensurate  [${\bf q}\ne(0\,,0)$] instability toward the $d$PI phase. 
The commensurate phase is sketched in \fig{fig1} (b) in real space. 
In momentum space it leads to Fermi surface deformations where 
the Fermi surface expands along the $k_y$ direction and shrinks 
along the $k_x$ direction [\fig{fig1} (e)], 
or vice versa if the bond along the $y$ direction would become stronger than the $x$ direction. 
While the commensurate $d$PI has been discussed since 2000,\cite{yamase00,metzner00} 
an incommensurate $d$PI starts to be discussed 
very recently.\cite{metlitski10,metlitski10a,holder12,husemann12} 

c) $V^a = (0,0,1,0,0,0)$ ($(0,0,0,1,0,0)$), which correspond to the freeze of 
the third (fourth) component and describe the instability toward 
the BOP$_x$ (BOP$_y$).\cite{cappelluti99,foussats04,morse91} 
The corresponding modulation vector turns out to be ${\bf q}=(\pi,\pi)$  or 
very close to it. The commensurate BOP$_x$, namely with ${\bf q}=(\pi,\pi)$, 
is sketched in \fig{fig1} (c) 
whereas the BOP$_y$ with ${\bf q}=(\pi,\pi)$ 
is obtained by rotating \fig{fig1} (c) by 90$^{\circ}$.  

d) $V^a = \frac{1}{\sqrt{2}} (0,0,1,1,0,0)$, which corresponds to the freeze of 
both third and fourth components simultaneously. 
The modulation vector is ${\bf q}=(\pi,\pi)$ or 
very close to it, as in the case of the BOP$_{x}$ and BOP$_{y}$. 
We refer to this instability as the BOP$_{xy}$. 
The commensurate BOP$_{xy}$ with ${\bf q}=(\pi,\pi)$ 
is sketched in \fig{fig1} (d). For simplicity, we also use the phrase BOP 
when we do not have to distinguish between BOP$_x$, BOP$_y$, and BOP$_{xy}$.


e) $V^a = (1,0,0,0,0,0)$, which corresponds to the freeze of 
charge fluctuations $\delta R$ and describes 
the instability toward PS for ${\bf q}=(0,0)$ and a charge-density-wave phase, 
including stripes, for a finite ${\bf q}$. 
A finite ${\bf q}$ instability, however, was not detected in the 
present study. 


In general, eigenvectors of $[D_{ab}({\bf q}, \rm{i}\omega_n)]^{-1}$ can have a non-zero 
value in each component. 	
However we checked that the inner product
between an eigenvector of $[D_{ab}({\bf q}, \rm{i}\omega_n)]^{-1}$ and $V^a$ 
becomes larger than $0.99$ at the corresponding critical point. 
 
\subsection{Effective susceptibilities}
While numerical results presented in this paper are 
computed from the full susceptibility \eq{dyson}, 
it is instructive to extract an effective susceptibility associated with 
each instability explained in the previous section 
by discarding the interactions with other modes contained in $D_{ab}({\bf q}, \rm{i}\omega_n)$. 

The usual charge-charge correlation function is written as 
\begin{eqnarray}\label{chi}
\chi^c_{ij}(\tau)=\frac{1}{N} \sum_{p,q} \langle T_\tau X_i^{pp}(\tau)
X_j^{qq}(0)\rangle \; .
\end{eqnarray}
Using the completeness condition [\eq{eq:v1}] and the relation between
$X^{00}_i$ and $\delta R_i$ [\eq{lambdai}], $\chi^c_{ij}$ can be written  
in Fourier space,\cite{foussats04}
\begin{eqnarray}\label{chi-c}
\chi^c({\bf q},{\rm i}\omega_n)=
-N {\left( \frac{\delta}{2} \right)}^2 D_{11}({\bf q},{\rm i}\omega_n)\,.
\end{eqnarray}
Thus the charge-charge correlation function is just the component $(1,1)$ of the $D_{ab}$.
Note that the factor $N$ in front of the right hand side of \eq{chi-c} 
shows that 
charge fluctuations are of $O(1)$ since $D_{ab} \propto 1/N$ as seen in \eq{D0inverse}. 

The susceptibility of the $d$CDW is obtained by 
focusing on the sector $a, b=5, 6$ of the matrix $D_{ab}^{-1}$. 
We obtain 
\begin{eqnarray}
\chi_{d{\rm CDW}}({\bf q},{\rm i}\omega_n)=
 [(8/J) \Delta^2-\Pi_{d{\rm CDW}}({\bf q},{\rm i}\omega_n)]^{-1}\,,
\label{chi-flux}
\end{eqnarray}
where $\Pi_{d{\rm CDW}}({\bf q},{\rm i}\omega_n)$ is the electronic polarizability of 
the $d$CDW and is given by 
\begin{eqnarray}
\Pi_{d{\rm CDW}}({\bf q}, {\rm i}\omega_n) = - \frac{1}{N_{s}}\;
\sum_{{\bf k}}\; \gamma_{d{\rm CDW}}^2({\bf q},{\bf k}) \frac{n_{F}(\epsilon_{{\bf k} + {\bf q}}) 
- n_{F}(\epsilon_{{\bf k}})} 
{\epsilon_{{\bf k} + {\bf q}} - \epsilon_{{\bf k}}-{\rm i} \omega_n}\,,
\label{PiCDW}
\end{eqnarray}
with a form factor 
$\gamma_{d{\rm CDW}}({\bf q},{\bf k})=2 \Delta [\sin(k_x+q_x/2)-\sin(k_y+q_y/2)]$. 


Similarly, the susceptibility of the $d$PI is obtained from the sector 
$a, b=3, 4$ of the matrix $D_{ab}^{-1}$:   
\begin{eqnarray}
\chi_{d\rm PI}({\bf q},{\rm i}\omega_n)= [(8/J) \Delta^2-\Pi_{d\rm PI}({\bf q},{\rm i}\omega_n)]^{-1}
\label{chi-dPI}
\end{eqnarray}
and the electronic polarizability  of the $d$PI reads 
\begin{eqnarray}
\label{Pi-dPI}
\Pi_{d\rm PI}({\bf q}, {\rm i}\omega_n) = - \frac{1}{N_{s}}\;
\sum_{{\bf k}}\; \gamma_{d\rm PI}^2({\bf q},{\bf k}) \frac{n_{F}(\epsilon_{{\bf k} + {\bf q}}) 
- n_{F}(\epsilon_{{\bf k}})} 
{\epsilon_{{\bf k} + {\bf q}} - \epsilon_{{\bf k}}-{\rm i}\omega_n} \,,
\end{eqnarray}
with a form factor 
$\gamma_{d\rm PI}({\bf q},{\bf k})=2 \Delta [\cos(k_x+q_x/2)-\cos(k_y+q_y/2)]$. 

For the case of the BOP$_x$ we focus on the sector $a=b=3$ and obtain 
\begin{eqnarray}
\chi_{{\rm BOP}_x}({\bf q},{\rm i}\omega_n)= 
[(4/J)\Delta^2 -\Pi_{{\rm BOP}_x}({\bf q},{\rm i}\omega_n)]^{-1}\,,
\label{chi-BOP}
\end{eqnarray}
where the electronic polarizability is given by 
\begin{eqnarray}\label{Pi-BOP}
\Pi_{{\rm BOP}_x}({\bf q}, {\rm i}\omega_n) = - \frac{1}{N_{s}}\;
\sum_{{\bf k}}\; 4 \Delta^2 \cos^2 (k_x+q_x/2) \frac{n_{F}(\epsilon_{{\bf k} + {\bf q}}) 
- n_{F}(\epsilon_{{\bf k}})} 
{\epsilon_{{\bf k} + {\bf q}} - \epsilon_{{\bf k}}-{\rm i}\omega_n}\,. 
\end{eqnarray}
For the case of $a=b=4$, i.e., BOP$_y$, 
the form factor in \eq{Pi-BOP} is replaced by $\cos^2( k_y+q_y/2)$.  
It is easily seen in \eq{Pi-BOP} that the BOP$_{x}$ and BOP$_{y}$ 
instabilities occur simultaneously, but with a different modulation vector: 
suppose ${\bf q}=(q_x,q_y)$ for the BOP$_{x}$, then ${\bf q}=(q_y,q_x)$ 
for the BOP$_{y}$. While we will not present results for the BOP$_{y}$, 
it should be understood that the instability of the  BOP$_{y}$ also exists. 

The susceptibility associated with BOP$_{xy}$ is given by the same equation 
as \eq{chi-dPI}, except that 
the form factor $\gamma_{d{\rm PI}}$  in \eq{Pi-dPI} is replaced by 
$\gamma_{{\rm BOP}_{xy}} ({\bf q},{\bf k})=2 \Delta [\cos(k_x+q_x/2)+\cos(k_y+q_y/2)]$. 

The form factor $\gamma_{d\rm CDW}({\bf q},{\bf k})$  [$\gamma_{d\rm PI}({\bf q},{\bf k})$] 
has a ${\bf k}$ dependence of $\cos k_x-\cos k_y$ at ${\bf q}=(\pi,\pi)$ [${\bf q}=(0,0)$], 
which indicates the $d$-wave character of the instability. 
Note that the $d$PI and $d$CDW belong to different 
eigenspace and are not connected with each other by changing the momentum ${\bf q}$. 

While the terminology of the $d$PI itself makes sense when a modulation vector 
is close to ${\bf q}=(0,0)$, we may consider formally a large ${\bf q}$ 
in Eqs.~(\ref{chi-dPI}) and (\ref{Pi-dPI}). 
The $d$PI is then connected with the BOP$_{x}$ and  BOP$_{y}$ 
when ${\bf q}$ is located along the direction of $(\pi,0)$-$(\pi,\pi)$ or 
$(0,\pi)$-$(\pi,\pi)$. Suppose ${\bf q}'=(\pi,q_y)$, we can easily find 
\begin{equation}
\Pi_{d\rm PI}({\bf q}', {\rm i}\omega_n) =
\Pi_{{\rm BOP}_x}({\bf q}', {\rm i}\omega_n) + 
\Pi_{{\rm BOP}_y}({\bf q}', {\rm i}\omega_n)\,,
\end{equation}
by noting that 
\begin{equation}
\frac{1}{N_{s}} \sum_{{\bf k}}\; \sin k_x \cos (k_y + q_y /2) \frac{n_{F}(\epsilon_{{\bf k} + {\bf q}'}) 
- n_{F}(\epsilon_{{\bf k}})} 
{\epsilon_{{\bf k} + {\bf q}'} - \epsilon_{{\bf k}}-{\rm i}\omega_n}=0\,.
\end{equation} 
We thus obtain 
\begin{equation}
\chi_{d\rm PI}^{-1}({\bf q}',{\rm i}\omega_n)= 
\chi_{{\rm BOP}_x}^{-1}({\bf q}',{\rm i}\omega_n) + 
\chi_{{\rm BOP}_y}^{-1}({\bf q}',{\rm i}\omega_n) \,.
\label{dPI-BOP}
\end{equation}
In particular, when ${\bf q}'$ is equal to ${\bf Q}\equiv (\pi,\pi)$, \eq{dPI-BOP} is reduced to 
\begin{equation}
\chi_{d{\rm PI}}({\bf Q}, {\rm i}\omega_n)=\frac{1}{2}\chi_{{\rm BOP}_x}({\bf Q}, {\rm i}\omega_n)\,,
\label{dPI-BOP-pi}
\end{equation} 
because 
$\chi_{{\rm BOP}_x}({\bf Q}, {\rm i}\omega_n)=\chi_{{\rm BOP}_y}({\bf Q}, {\rm i}\omega_n)$. 
Similarly, we can also obtain 
\begin{equation}
\chi_{d{\rm PI}}({\bf Q}, {\rm i}\omega_n)=\chi_{{\rm BOP}_{xy}}({\bf Q}, {\rm i}\omega_n)\,.
\label{dPI-BOPxy}
\end{equation}
Hence when the static BOP susceptibility diverges at ${\bf q}=(\pi,\pi)$, 
the $d$PI susceptibility also diverges simultaneously at the same momentum 
unless it already diverges at a different momentum.  
In fact, the $d$PI with ${\bf q}=(\pi,\pi)$ is equivalent to the 
BOP$_{xy}$ and is interpreted as superposition of the 
BOP$_x$ and BOP$_y$ as seen in \fig{fig1} (d). 


\section{Results}
We choose the parameters, $J/t=0.3$ and $V/t=0.5$. We set $t=1$, and 
all quantities with dimension of energy are in units of $t$. 
Irrespective of the presence of the $V$-term, 
our theory catches intrinsic charge instabilities in the 2D $t$-$J$ model 
such as $d$CDW, $d$PI, and BOP, which are driven by the $J$-term.
We compute the full susceptibility \eq{dyson} for various choices of $t'$ 
by assuming the paramagnetic phase 
and determine possible charge instabilities in the plane of 
hole density $\delta$ and temperature $T$. 
At half-filling an analytical solution is  obtained and 
the $d$CDW, dPI, and BOP have the same onset temperature, 
$T_c=J/8=0.0375t$, at which the static field $\Delta$ [\eq{staticDelta}] also sets in. 
Away from half-filling ($\delta > 0.004$) our computation is fully numerical. 
Since we determine critical lines by studying the susceptibility, 
the transition is continuous. In other words, a possibility of a first order 
transition is not considered in the present analysis. 

\subsection{Results for \boldmath{$t'=0$}}
Figure~\ref{fig2} shows the phase diagram for $t'=0$.
As mentioned in Sec.~II.B, five different types of charge instability are 
found: $d$CDW, $d$PI, BOP$_x$, 
BOP$_{xy}$, and PS. 
The instability toward the commensurate $d$CDW, namely with ${\bf q}=(\pi,\pi)$, 
occurs in a wide doping region. The transition temperature decreases 
gradually with increasing hole density and exhibits reentrant behavior at low $T$ 
in the region $0.12 \lesssim \delta \lesssim 0.14$. 
However, near $\delta \approx 0.14$  an incommensurate [${\bf q}\ne(\pi,\pi)$] $d$CDW instability 
occurs below $T\approx 0.015$ and its critical doping rate is higher than that of 
the commensurate $d$CDW. Hence the resulting critical line of the $d$CDW 
follows the outer line, i.e., the thin line at low $T$ and the thick line for high $T$.

\begin{figure}[htp]
 \centering
 \includegraphics[width=8cm]{figure2.eps}
 \caption{(Color online) Critical temperature versus doping rate 
for $d$CDW, $d$PI, BOP$_x$, BOP$_{xy}$, and PS for $t'=0$ and $J=0.3$. 
Thick (thin) lines describe commensurate (incommensurate) transitions. 
The critical line for PS is shown in a larger 
scale of $T$ in the inset.}
\label{fig2}
\end{figure}



For $t'=0$ 
the commensurate $d$PI with ${\bf q}=(0,0)$ has the same onset temperature as 
the commensurate BOP$_x$ and BOP$_{xy}$, namely with ${\bf q}=(\pi,\pi)$; 
this reason will be explained in the last paragraph in the present section. 
The transition lines exhibit reentrant behavior at low $T$.   
However, an incommensurate $d$PI  
emerges at low $T$
and preempts the reentrant line, extending the region of $d$PI. 
BOP$_x$ and BOP$_{xy}$ also exhibit an incommensurate instability 
at low $T$ and preempt their commensurate instabilities. 
Furthermore the degeneracy between BOP$_x$ and BOP$_{xy}$ is lifted 
via an incommensurate transition. 
While the critical doping rate for BOP becomes  higher than for $d$PI at low $T$, 
this result occurs only for a small $t'$ and, as will be shown below, 
the opposite occurs in the presence of a realistic $t'$ for cuprates.  

The system also exhibits PS  at low doping. 
The inset of Fig. 2 shows the PS line in a larger temperature scale. 
We see that PS appears with decreasing temperature, 
but with further cooling down it goes back to the 
paramagnetic phase. As a result, PS occurs only 
in an intermediate temperature region. This peculiar reentrant behavior 
was also found in Ref.~\onlinecite{koch04} in the Hubbard model. 
The region of PS shrinks with increasing $V$ and also by introducing $t' (<0)$, 
as will be discussed in the subsection D.  


The phase diagram in \fig{fig2} should not 
be interpreted in such a way that the $d$CDW is unstable against 
the $d$PI or BOP at low $T$ or low $\delta$, 
because we perform a stability analysis 
in the paramagnetic phase. 
Rather, \fig{fig2} indicates a hierarchy for different charge instabilities, i.e., 
the outer the critical line is, the stronger the tendency toward 
the corresponding instability is. 


It requires highly accurate numerics to determine precisely a modulation vector ${\bf q}$ 
of each order along its outer critical line because of a rather flat structure of the susceptibility 
in momentum space, especially for the $d$PI. 
Therefore considering our achieved numerical accuracy we present in \fig{fig3}   
modulation vectors of each instability, for which the absolute value of the 
corresponding eigenvalue of 
$D_{ab}^{-1}({\bf q},0)$ becomes less than $10^{-4}t$ on its outer critical line. 
The width of such a ${\bf q}$ region, at a fixed temperature, implies 
how sharp the susceptibility is in momentum space. 
Since modulation vectors of each instability are computed along its outer critical line, 
each critical temperature shown in \fig{fig3} corresponds to a certain critical doping rate, 
which can be read off from \fig{fig2}. 
Although \fig{fig3} is presented only along the axis $(\pi,0)$-$(\pi,\pi)$-$(0,0)$-$(\pi,0)$, 
we scanned the  whole ${\bf q}$ region of the Brillouin zone and checked numerically 
that instabilities indeed occur along that axis. 


\vspace{5mm}
\begin{figure}[htp]
 \centering
 \includegraphics[width=12cm]{figure3.eps}
\caption{(Color online) 
Modulation vectors for $d$CDW (left panel), $d$PI (middle panel), and 
BOP$_x$ (right panel) along the corresponding outer critical line 
in \fig{fig2}. For each critical temperature, the critical doping rate can be 
read off from \fig{fig2}.  Because of symmetry,  the results along the 
direction of $(0,\pi)$-$(\pi,\pi)$ are the same as those along  
the $(\pi,0)$-$(\pi,\pi)$ direction for the $d$CDW and $d$PI. 
}
\label{fig3}
\end{figure}



In \fig{fig3} (a) we show the result for the $d$CDW along its outer critical line in \fig{fig2}. 
At high critical temperature (i.e., low critical doping rate) the instability occurs at ${\bf q} = (\pi,\pi)$,
and with lowering temperature the modulation vector shifts 
from $(\pi,\pi)$ and becomes incommensurate. 
In \fig{fig3} (c) we plot the corresponding modulation vector of the BOP$_x$. 
At low $T$, the modulation vector ${\bf q}$ slightly shifts 
from $(\pi,\pi)$ and the BOP$_x$ becomes incommensurate, as in the case of the $d$CDW.  
At high $T$, ${\bf q}$ is located at $(\pi,\pi)$, but in contrast to the case of the 
$d$CDW, the ${\bf q}$ region is not extended on the side of the direction of 
$(\pi,\pi)$-$(0,0)$. This is because the eigenvector of the BOP$_x$ does not exist there, 
instead, the eigenvector of the full susceptibility [\eq{dyson}] 
changes to that corresponding to the $d$PI. 
A modulation vector of the BOP$_{xy}$ appears only along the axis 
$(\pi,\pi)$-$(0,0)$. Its $T_c$ dependence is very similar to that of BOP$_x$, but 
the labels $(\pi,0)$ and $(0,0)$ in \fig{fig3} (c) 
should be replaced by $(0,0)$ and $(\pi,0)$, respectively. 


The corresponding result for the $d$PI is shown in \fig{fig3} (b), which looks 
very different from those for the $d$CDW and BOP. 
In fact, for $t'=0$, the static electronic polarizability of the $d$PI 
has a special feature, which was already noted in Ref.~\onlinecite{morse91} 
in a different context. 
To see this we rewrite \eq{Pi-dPI} in a different form, 
\begin{equation}
\Pi_{d\rm PI}({\bf q}, 0) = - \frac{1}{N_{s}}\;
\sum_{{\bf k}}\;   n_{F}(\epsilon_{{\bf k}})\left[
\frac{\gamma_{d\rm PI}^2({\bf q},{\bf k})}{\epsilon_{{\bf k}}-\epsilon_{{\bf k}+{\bf q}}} 
+\frac{\gamma_{d\rm PI}^2(-{\bf q},{\bf k})}{\epsilon_{{\bf k}}-\epsilon_{{\bf k}-{\bf q}}} 
\right]\,.
\end{equation}
When ${\bf q}$ lies along the diagonal direction ${\bf q} \parallel (q,q)$, 
we find after some algebra 
\begin{equation}
\Pi_{d\rm PI}({\bf q}, 0) = \frac{8 \Delta^2}{(t\delta+2\Delta) N_{s}}\;
\sum_{{\bf k}}\;   n_{F}(\epsilon_{{\bf k}})\cos\frac{k_x+k_y}{2}
\tan\frac{k_x-k_y}{2}\sin\frac{k_x-k_y}{2}\,, 
\label{d-bubble}
\end{equation}
that is, the static $d$PI susceptibility [\eq{chi-dPI}] 
does not depend on ${\bf q}$ for any momentum along the diagonal direction. 
This result holds for any carrier density and any temperature. 
Therefore, if the $d$PI takes place for a vector ${\bf q}$ in the diagonal direction, 
the susceptibility diverges simultaneously at all ${\bf q}$ along the diagonal direction. 
The full susceptibility [\eq{dyson}] actually shows that feature in \fig{fig3} (b). 
Furthermore,  this flat feature of the susceptibility extends more away from the diagonal direction. 
The ${\bf q}$ region along $(0,0)$-$(\pi,0)$ shrinks at $T_c \approx 0.024$ in \fig{fig3} (b), 
which results from the proximity to PS, as will be discussed in the subsection~C. 
While the susceptibility is always flat along the diagonal direction of ${\bf q}$, 
the susceptibility shows a peak at a modulation vector along 
$(0,0)$-$(\pi,0)$ at low $T$. 
The ${\bf q}$ region has a sharp boundary at $(\pi,\pi)$ in \fig{fig3} (b) and 
the $d$PI does not have any possible modulation vector 
along the $(\pi,0)$-$(\pi,\pi)$ direction. 
This is because 
the eigenvector corresponding to the $d$PI is not realized along 
$(\pi,0)$-$(\pi,\pi)$, instead,  the BOP$_x$ eigenvector appears there. 
This property may be understood also 
in terms of the effective susceptibilities. 
Equation~(\ref{dPI-BOP}) indicates that if $\chi_{d\rm PI}^{-1}({\bf q}',0)$ becomes 
zero, either $\chi_{{\rm BOP}_x}^{-1}({\bf q}',0)$ or $\chi_{{\rm BOP}_y}^{-1}({\bf q}',0)$ 
should be already negative,  
since in general $\chi_{{\rm BOP}_x}$ is not equal to $\chi_{{\rm BOP}_y}$ 
for a momentum along $(\pi,0)$-$(\pi,\pi)$, except for 
${\bf q}'=(\pi,\pi)$ where both $\chi_{{\rm BOP}_x}^{-1}({\bf q}',0)$ and 
$\chi_{{\rm BOP}_y}^{-1}({\bf q}',0)$ can become zero simultaneously.   
Therefore a possible instability of the $d$PI along $(\pi,0)$-$(\pi,\pi)$ 
is replaced by the BOP$_x$. 


The ${\bf q}$-independence of $\Pi_{d{\rm PI}}$ along the diagonal direction 
leads to another special feature. As we mentioned at the end of Sec.~II, 
the onset temperature of the BOP with ${\bf q}=(\pi,\pi)$ is the same as that of the $d$PI 
with ${\bf q}=(\pi,\pi)$ [Eqs.~(\ref{dPI-BOP-pi}) and (\ref{dPI-BOPxy})]. 
Therefore the onset temperature of 
the commensurate BOP becomes the same as that of the 
$d$PI with ${\bf q}={\bf 0}$, as shown in \fig{fig2}. 


\subsection{Results for finite \boldmath{$t'$}}
The degeneracy between the $d$PI and BOP seen in \fig{fig2} 
is lifted by introducing $t'$. 
The upper and lower panels in \fig{fig4} show 
the results for $t'=-0.20$ and $-0.30$, respectively. 
While the BOP instability is always restricted to a lower doping region,  
the $d$PI becomes favorable in a wider doping region. 
Near half-filling the $d$PI and BOP are still almost degenerate because,    
as seen from \eq{Ek}, the hopping integral $t'$ is renormalized to be 
$t' r_0 \propto t' \delta$ and becomes irrelevant close to half-filling. 
The BOP$_x$ and BOP$_{xy}$ are always degenerate as far as 
they exhibit a commensurate instability. 
Their degeneracy is lifted when their modulation 
vector becomes incommensurate at low temperature. 

\begin{figure}[tp]
 \centering
 \includegraphics[width=8cm]{figure4.eps}
 \caption{(Color online) Critical temperature and doping rate 
for $t'=-0.20$ (upper panel) and $t'=-0.30$ (lower panel). 
The notation is the same as \fig{fig2}. 
}
 \label{fig4}
\end{figure}


As shown in \fig{fig4}, the doping region of the commensurate $d$CDW instability 
is extended by the presence of $t'$ and 
an incommensurate $d$CDW becomes dominant 
at high $\delta$ and low critical temperature. 
On the other hand, PS is suppressed by introducing $t'$. 
The critical line for PS bends back to zero doping for high $T$ (not shown) in a way similar 
to the case for $t'=0$ (inset of \fig{fig2}). 

The modulation vector of each instability is shown in the upper 
and middle row in \fig{fig5} for $t'=-0.20$ and $-0.30$, respectively, 
along the corresponding outer critical line in \fig{fig4}. 
Both $d$CDW [Figs.~\ref{fig5} (a) and (d)] and BOP$_x$ [Figs.~\ref{fig5} (c) and (f)] 
show an instability 
at ${\bf q}=(\pi,\pi)$ for high critical temperature, and becomes incommensurate for low critical temperature. 
These features are the same as those seen in Figs.~\ref{fig3} (a) and (c). 
Results of BOP$_{xy}$ are essentially the same as those of BOP$_{x}$, 
but the labels $(\pi,0)$ and $(0,0)$ should be interchanged in Figs.~\ref{fig5} (c) and (f). 


\begin{figure}[tp]
 \centering
 \includegraphics[width=12cm]{figure5.eps}
\caption{(Color online) Modulation vectors for $d$CDW (left panels), $d$PI (middle panels), 
and BOP$_x$ (right panels) for $t'=-0.20$ (upper row), $t'=-0.30$ (middle row),
and $t'=-0.40$ (lower row) along the outer critical line of the corresponding order 
for each $t'$ in \fig{fig4} and \fig{fig6} (middle).  
}
 \label{fig5}
\end{figure}


Figs.~\ref{fig5} (b) and (e) show results for the $d$PI. 
They are very different from \fig{fig3} (b), 
except for a region of high critical temperature near $T_c\approx 0.03$, 
namely the doping region $0< \delta_c \lesssim 0.02,$ where 
the effect of $t'$ becomes irrelevant. 
For $t'=-0.20$ [\fig{fig5} (b)], as the critical temperature decreases, 
the ${\bf q}$ region shrinks around ${\bf q}=(0,0)$ and 
the $d$PI tends to become commensurate. 
Close to zero temperature, however, 
a tendency toward an incommensurate $d$PI appears 
in the $(0,0)$-$(\pi,0)$ direction. 
This incommensurate feature is also seen more clearly for $t'=0$ at low $T$ [\fig{fig3} (b)], 
and disappears quickly with increasing $t'$. It becomes 
nearly invisible for $t'=-0.30$ in our temperature scale. 
For $t'=-0.30$,  
the ${\bf q}$ region shrinks first around ${\bf q}=(0,0)$ with 
decreasing the critical temperature. 
In the intermediate temperature range, $0.021 \gtrsim T_c \gtrsim 0.008$, 
the modulation vector becomes slightly incommensurate along the diagonal.
This deviation from the commensurate vector is also barely visible 
in the bottom panel of \fig{fig4},  where the incommensurate 
$d$PI line separates very slightly from the commensurate $d$PI in the corresponding 
temperature region. For $T_c \lesssim 0.008$, the $d$PI becomes fully commensurate. 


\begin{figure}[tp]
 \centering
 \includegraphics[width=12cm]{figure6.eps}
 \caption{(Color online) The same plot as \fig{fig4}, but for different choices of $t'$: 
 $t'=-0.35$ (top panel), $t'=-0.40$ (middle panel), and $t'=-0.45$ (bottom panel). 
 In the latter two cases the phase diagram in a high doping region is 
 magnified in the inset.   
          }
 \label{fig6}
\end{figure}

The phase diagram close to half-filling 
does not depend essentially on a choice of $t'$.  
In fact, the critical lines for the 
BOP and PS do not change much even for a further larger $t'$. 
However, we find that tendencies toward $d$CDW and $d$PI have strong $t'$ 
dependence. 
In \fig{fig6} we present the phase diagram for $t'=-0.35$, $-0.40$, and $-0.45$. 
First we focus on the $d$CDW for $t'=-0.35$ and $-0.40$. 
The instability extends to a higher doping region with increasing $t'$ 
and the commensurate $d$CDW tends to become more favorable even at low $T$.  
The doping region of the $d$CDW, however,  
starts to decrease quickly for $|t'| > 0.41$, as can be seen in the result for $t'=-0.45$. 


In contrast to the $d$CDW, the outer critical line for the $d$PI 
extends to higher doping with increasing $t'$. 
For $t'=-0.35$, an incommensurate $d$PI becomes dominant 
in a wide doping region ($0.05 \lesssim \delta \lesssim 0.25$). 
While the commensurate $d$PI is realized for 
$\delta \lesssim 0.16$ and $\delta \gtrsim  0.19$, 
it does not occurs between these two doping region. 
This feature is more evident for $t'=-0.40$ [\fig{fig6} (middle panel)]. 
In a wide doping region ($0.10 \lesssim \delta \lesssim 0.30$) 
the commensurate $d$PI does not occur 
and only an incommensurate $d$PI is possible. 
However in a high doping region ($0.30 \lesssim \delta \lesssim 0.35$), 
the commensurate $d$PI shows up again; 
see the inset in  \fig{fig6} (middle panel). 
These peculiar features of the $d$PI are due to the presence of the van Hove singularity around 
$\delta = 0.33$, where the $d$-wave weighted density of states,\cite{yamase05}  
which is defined by $\lim_{{\bf q} \rightarrow {\bf 0}} \Pi_{d\rm PI}({\bf q},0)$ in \eq{Pi-dPI},  
is enhanced, favoring the instability at ${\bf q}=(0,0)$. 
In an intermediate doping region, 
the $d$-wave density of states 
is suppressed and the $d$PI with a finite ${\bf q}$ 
becomes more favorable. 
Closer to half-filling, however, the $d$-wave density of states is again 
enhanced because of narrowing the band width upon approaching half-filling, 
leading to a recovery of the commensurate 
$d$PI for $\delta \lesssim 0.10$. 
For $t'=-0.45$ these features are more emphasized. 
Figure~\ref{fig6} (bottom panel) shows that 
the commensurate $d$PI occurs both, close to half-filling and 
around van Hove filling ($\delta=0.42$), and these two regions are connected 
by the critical line of an incommensurate $d$PI. 
The commensurate $d$PI close to half-filling and 
in a high doping region was also found in the slave-boson mean-field theory.\cite{yamase00} 

In \fig{fig5} we compare the modulation vector of each instability for 
$t'=-0.20$, $-0.30$, and $-0.40$. 
With increasing $t'$, the $d$CDW tends to be more commensurate 
even at low critical temperature. 
For the $d$PI, on the other hand, a large $t'$ tends to favor 
an incommensurate modulation along 
the diagonal direction of the Brillouin zone in an intermediate range of a critical 
temperature, and the commensurate $d$PI can be realized only at high and low critical temperature, 
corresponding to a doping region close to half-filling and around van Hove filling, 
respectively. The modulation vectors of the BOP$_x$ (and also BOP$_{xy}$) 
do not depend much on a choice of $t'$. 


\subsection{Mutual interaction among different modes}
Because of the renormalization of the bosonic propagators due to the coupling 
to electronic bubbles [\eq{dyson}], one would expect in general some coupling 
among different modes. 
This effect actually appears for modulation vectors of the $d$PI. 
There is a dip on the side of the region of 
$(0,0)$-$(\pi,0)$ at $T_c\approx 0.024$ in \fig{fig3} (b), and 
$T_c \approx 0.025$ and $0.026$ in Figs.~\ref{fig5} (b) and (e), respectively. 
This dip occurs near the temperature where the critical lines of the $d$PI and PS 
cross each other (see Figs.~\ref{fig2} and \ref{fig4}). 
Moreover, we checked that the dip in question does not appear in results obtained 
from the effective susceptibility of the $d$PI [\eq{chi-dPI}]. 
Therefore a tendency toward PS plays a role in the suppression of the 
incommensurate feature of the $d$PI along $(0,0)$-$(\pi,0)$. 

Except for the above feature, we checked that our results (Figs.~\ref{fig2}-\ref{fig6})  
are nearly the same as those determined by the effective susceptibilities 
[Eqs.~(\ref{chi-c}), (\ref{chi-flux}), (\ref{chi-dPI}), and (\ref{chi-BOP})]. 
In this sense, the coupling among different bosonic fluctuations is rather weak 
at least in leading order.



\subsection{Effect of the \boldmath{$V$}-term and stability of phase separation} 
We checked that the results for $d$CDW, $d$PI, and BOP are almost intact 
for different choices of $V$ $(\geq 0)$ 
and  that an additional instability such as the usual checkerboard charge density wave 
does not occur at least for $V \leq 1$ for 
any doping rate.\cite{foussats04,misc-foussats04} 
Furthermore, the reentrant critical line of PS (inset of \fig{fig2}) is a robust feature. 
However, it is a subtle issue whether PS actually occurs at $T=0$. 
The result depends on choices of $V$, $t'$, and $J$. 
For $t'=0$ and $J=0.3$, we found that PS at $T=0$ occurs for $\delta \lesssim 0.08$ $(0.025)$ 
at $V=0$ $(0.1)$ and vanishes already for $V\gtrsim 0.2$. 
When $t'$ is introduced, PS is strongly suppressed even at $V=0$, 
for example, it occurs at $T=0$ only for $\delta \lesssim 0.01$  for $t'=-0.35$. 
A smaller $J$ also suppresses PS, and for the special case of $J=0$ 
no PS is observed at $T=0$ for any $V$ $(\geq 0)$ and $t'$ $(\leq 0)$.  


\section {Conclusion} 
Applying a large-$N$ expansion formulated in the path integral representation 
of the $t$-$J$ model, we have analyzed all possible 
charge instabilities of the paramagnetic phase, and have elucidated  
the phase diagram in the doping and temperature plane for 
a sequence of $t'$.
We have found that $d$CDW, $d$PI, BOP, and PS are the most important 
charge instabilities. 
The first two instabilities appear in a wide doping region. 
The $d$CDW usually becomes the leading instability and the $d$PI 
occurs as a next leading one with a strong 
tendency to become incommensurate. 
In the presence of a large $t'$, however, 
we have found that the $d$PI becomes the leading instability 
in a high doping region.  
Considering the high complexity of the $t$-$J$ model, it 
is beyond the scope of the present study to address 
which charge instability would become ultimately the leading one. 
Rather, the present stability analysis on charge instabilities of the paramagnetic state 
was motivated by active studies on the PG in cuprates in terms of 
various charge fluctuations, and we have provided a microscopic basis of possible charge fluctuations in doped Mott insulators. 
We first compare our results with literature and then discuss implications for 
cuprate superconductors. 

Taking into account a number of papers about charge stripes 
in cuprates,\cite{kivelson03,vojta09} it may be surprising that 
we do not detect  stripe tendencies in our formalism, 
which  exclusively favors charge instabilities. If the $t$-$J$ model would exhibit  
a tendency to charge stripes order, it might be a consequence of 
a coupling with incommensurate magnetic modes, 
which could appear in the next-to-leading order [$O(1/N)$] in the present scheme. 
A charge-stripe solution was actually obtained in the $t$-$J$ model 
with $t'=0$ in the DMRG study,\cite{white98,white00,corboz11} which 
however contradicts other studies.\cite{hellberg99,hu12} 
In the presence of $t'$, on the other hand, 
most of numerical studies in the $t$-$J$ model 
reported that the charge stripe solution becomes unstable.\cite{white99,tohyama99} 
Our results, therefore, agree to major literature.  

Interestingly, the BOP$_{x (y)}$ shares the same feature as stripe order 
from a symmetry point of view. 
When ${\bf q}$ shifts away from $(\pi,\pi)$, 
${\bf q}$ is located only along the $q_{y(x)}$ direction and thus 
the BOP$_{x(y)}$ breaks both orientational and translational symmetry of the lattice. 
In fact such an incommensurate BOP$_{x(y)}$ instability is found to occur 
up to $\delta \sim 0.10$;  see Figs.~\ref{fig2},~\ref{fig4},~and \ref{fig6}. 
While the BOP has not been discussed much so far, 
the BOP was also obtained in other studies in 
the $t$-$J$ model.\cite{cappelluti99,foussats04,morse91} 

As discussed in Sec.~III.D, 
PS at $T=0$ strongly depends on a choice of $V$, $t'$, and $J$.  
However it is a robust feature that PS occurs in an intermediate temperature 
region.  This property for a finite $T$ has not been discussed so far 
except for Ref.~\onlinecite{koch04} in the Hubbard model, probably 
because various numerical simulations are 
usually coded only at zero temperature. 
Interestingly, the reentrant critical line of PS (see the inset of \fig{fig2}) 
was interpreted as a source to generate a strong forward scattering channel of 
the electron-phonon vertex which emerges only at finite $T$.\cite{koch04} 

While the commensurate $d$PI (${\bf q}={\bf 0}$) was already found  in the 
$t$-$J$\cite{yamase00} and Hubbard\cite{metzner00} models, 
an incommensurate $d$PI (${\bf q} \ne {\bf 0}$) started to be discussed 
quite recently.\cite{metlitski10,metlitski10a,holder12,husemann12} 
We have obtained that the static $d$-wave electronic polarizability [\eq{d-bubble}] 
does not depend on ${\bf q}$ for any momentum along the Brillouin zone diagonal, 
which holds for any temperature and any electron density as long as $t'$ is zero. 
In our model, this feature remains even for a finite $t'$ near half-filling 
since $t'$ is effectively reduced by a factor of $\delta$. 
Our result shown in the inset of \fig{fig6} (bottom panel) may be best 
compared with existing results,  since they were obtained 
in a weak coupling analysis.\cite{metlitski10,metlitski10a,holder12,husemann12} At hole density below van Hove filling ($\delta \lesssim 0.42$), 
the leading instability is an incommensurate $d$PI and its modulation vector 
is located along the $(0,0)$-$(\pi,\pi)$ direction, which 
agrees with literature.\cite{metlitski10,metlitski10a,holder12,husemann12} 
At hole density above van Hove filling we have obtained the $d$PI with ${\bf q}={\bf 0}$,  
while Ref.~\onlinecite{holder12} showed that the static electronic polarizability 
of the $d$PI has a peak 
along the $(0,0)$-$(\pi,0)$ direction at least at $T=0$. 
In our case, an incommensurate peak along the 
$(0,0)$-$(\pi,0)$ direction indeed develops as seen in Figs.~\ref{fig3}~(b), 
\ref{fig5}~(b) for a small $t'$, but it develops 
below extremely low temperature for a large $t'$. 
This effect is not visible in \fig{fig6} (bottom). 

In the $1/N$ expansion, the $d$CDW is the leading instability in most of cases, 
in agreement with previous studies.\cite{morse91,cappelluti99,foussats04} 
We have found that close to the $d$CDW, the $d$PI also exists in a wide doping region.  
Therefore fluctuations associated with 
both $d$CDW and $d$PI are expected to be important for temperatures 
above the onset of the $d$CDW. 
The mutual interaction between $d$CDW and $d$PI  
seems rather weak at least in leading order 
because both, the full calculation [\eq{dyson}] and effective 
susceptibilities [Eqs.~(\ref{chi-c}), (\ref{chi-flux}), (\ref{chi-dPI}), and (\ref{chi-BOP})], 
give nearly the same results. 
Given that the presence of $t'$ is usually assumed for cuprates 
and our critical lines of $d$CDW and $d$PI have the same 
doping dependence of the PG temperature, furthermore with a comparable temperature 
if $t \approx 500$ meV, 
it is interesting to study each fluctuation effect on the 
electronic property. In fact,  
existing work already showed interesting results, but 
with some open questions. 
The phenomenological study assuming the $d$CDW showed that essential 
features of the PG are well captured.\cite{chakravarty01} 
However fluctuation effects were not considered in Ref.~\onlinecite{chakravarty01}. 
In a perturbative calculation of the electron self-energy due to 
a coupling to $d$CDW fluctuations in the $t$-$J$ model, 
a pseudogap is obtained in the electronic spectral function, 
which shares many important features with experimental 
data.\cite{greco09,bejas11,greco11}  
Moreover the same analysis was also applied to the 
explanation of anisotropic scattering rate of quasi-particles\cite{buzon10} 
observed in angle-dependent magnetoresistance experiments.\cite{abdel06}  
However, no calculation was performed beyond the perturbative analysis. 
On the other hand, in a perturbative calculation for $d$PI fluctuations centered 
around ${\bf q}=(0,0)$,  
a splitting of the spectral function near the Fermi energy was obtained, reminiscent 
of a pseudogap.\cite{yamase12} 
Going beyond the perturbation theory and summing up all diagrams 
in the Gaussian fluctuation regime, however, instead of a splitting,  
the spectral function exhibits a broad single peak 
centered at the Fermi energy with a strong ${\bf k}$ dependence of 
$d$-wave symmetry.\cite{yamase12} The spectrum in the Ginzburg 
region\cite{ginzburg60} is an open question. 
Furthermore, the  role of incommensurate $d$PI fluctuations on a pseudogap 
phenomenon remains elusive. 

The $d$PI couples directly with $xy$ anisotropy such as due to a lattice structure 
and an external strain. While the $d$PI changes to a crossover phenomenon 
in such a case, the anisotropy can be strongly enhanced by the underlying $d$PI 
fluctuations as already discussed.\cite{yamase00,okamoto10} 
The same idea is also discussed 
for iron-pnictide superconductors near the structural phase transition 
from the tetragonal to orthorhombic phase.\cite{fisher11}  
Given that lattice anisotropy 
frequently exists in cuprates, 
the relevance of the $d$PI channel in the $t$-$J$ model suggests 
important implications for understanding cuprate superconductors,  
not only for Y-based\cite{hinkov07,hinkov08,daou10,yamase06,yamase09,hackl09} 
but also for La-based\cite{yamase00,yamase01,yamase02,yamase07} compounds. 

For a large $t'$, the commensurate $d$PI appears in a heavily overdoped region 
around van Hove filling ($\delta \approx 0.3-0.45$ in \fig{fig6}). 
While our critical line exhibits reentrant behavior at low $T$,  
the canonical model for the $d$PI\cite{khavkine04,yamase05} 
suggests that the reentrant behavior is preempted 
by a first order transition as a function of the chemical potential,  
or equivalently a phase separation as a function of doping, as far as 
the $d$PI occurs at ${\bf q}={\bf 0}$.  
It is known that Sr$_3$Ru$_2$O$_7$ exhibits the $d$PI in a magnetic 
field.\cite{grigera04,borzi07,rost09} 
In addition, a highly overdoped region in cuprates, 
where no superconducting and antiferromagnetic instabilities are expected, 
may also provide an opportunity to explore the $d$PI. 



\begin{acknowledgments}
The authors thank C. Gazza, T. Holder, C. Husemann, 
W. Metzner, S. Sorella, and R. Zeyher for valuable discussions. 
M. B. thanks Consejo Nacional de Investigaciones 
Cient\'{i}ficas y T\'ecnicas (CONICET) for financial support. 
A. G. thanks the National Institute for Materials Science (NIMS), 
where this work was initiated, and the Max-Planck-Institute for hospitality.
H. Y. was supported by a Grant-in-Aid for Scientific Research from Monkasho 
and the Alexander von Humboldt Foundation. 
\end{acknowledgments}


