text:
# Bond-order-wave phase and quantum phase transitions in the one-dimensional extended Hubbard model

## I Introduction

The one-dimensional (1D) extended Hubbard model has been extensively studied in recent years, both as an important theoretical test-bed for studying novel concepts in 1D (e.g., spin-charge separation) and methods, (e.g., quantum Monte Carlo, exact diagonalization, and the density matrix renormalization group) and as a useful model for several classes of quasi 1D materials including copper-oxide materials related to the high-Tc cuprate superconductors, conducting polymers  and organic charge-transfer salts. General 1D extended Hubbard models differ from the standard Hubbard model, which includes only an on-site electron-electron interaction ($`U`$), by the addition of longer-range interactions which are necessary to explain several experimentally observed effects in real materials, e.g., excitons in conducting polymers. The simplest extended Hubbard model (henceforth, EHM), on which we focus in this article, consists of adding a nearest neighbor interaction $`V`$. If the interaction parameters are assumed to arise solely from Coulomb interactions, both $`U`$ and $`V`$ are repulsive (positive), and $`U>V`$. However, viewed as phenomenological parameters incorporating the effects of additional (e.g. electron-phonon) interactions, the ranges of these parameters can be much broader, including $`U,V<0`$. The Hamiltonian is

$`H`$ $`=`$ $`t{\displaystyle \underset{i,\sigma }{}}(c_{i+1,\sigma }^{}c_{i,\sigma }+h.c.)`$ (3)
$`+U{\displaystyle \underset{i}{}}(n_{i,}\frac{1}{2})(n_{i,}\frac{1}{2})`$
$`+V{\displaystyle \underset{i}{}}(n_{i+1}1)(n_i1)\mu {\displaystyle \underset{i}{}}n_i,`$
where $`c_{i,\sigma }^{}(c_{i,\sigma })`$ creates (annihilates) an electron with spin $`\sigma `$ at site $`i`$, $`t`$ is the hopping integral between adjacent sites and $`\mu `$ is the chemical potential. Henceforth we set $`t=1`$ and express the interaction parameters $`U`$ and $`V`$ in units of $`t`$.

The ground state phase diagram of the EHM at half filling $`(\mu =0)`$ has been extensively studied using both analytical and numerical methods. Despite the apparent simplicity of the model, the phase diagram shows surprisingly rich structure. In the limit $`V`$=0 (the standard Hubbard model), the Hamiltonian (3) can be diagonalized exactly using the generalized Bethe Ansatz. For $`V0`$, the model has been studied using perturbative methods and numerical simulations. Broadly, the phase diagram consists of insulating phases with dominant charge-density-wave (CDW) and spin-density-wave (SDW) characters and metallic phases where singlet and triplet superconducting correlations dominate. In the physically relevant region for ”Coulomb-only” parameters ($`U,V>0`$), the system is in a CDW phase for large $`V/U`$ and in an state with dominant SDW fluctuations for small $`V/U`$. The CDW phase has broken discrete symmetry characterized predominantly by alternating doubly occupied and empty sites and exhibits long-range order. The SDW phase, on the other hand, has continuous symmetry and hence cannot exhibit long-range order in 1D (by the Mermin-Wagner theorem). Instead, it is a critical state characterized by the slow (algebraic) decay of the staggered spin-spin correlation function. Indeed, in the limit $`U1`$, $`UV`$, the model reduces to an effective Heisenberg model with $`J1/(UV)`$. For small $`U`$ and $`V`$ $`(U,V1)`$, the boundary between the CDW and the SDW phases was predicted to be at $`U=2V`$ using weak coupling renormalization group techniques (”g-ology”). Strong coupling calculations using second-order perturbation theory also gave the same phase boundary $`(U=2V)`$ between the CDW and the SDW phases for large $`U`$ and $`V`$ $`(U,V1)`$. For intermediate values of the parameters, the phase boundary was found to be shifted slightly away from the $`U=2V`$ line such that the SDW phase is enhanced, as shown by quantum Monte Carlo simulations  as well as strong coupling calculations using perturbation theory up to the fourth order. Moreover, the nature of the transition is quite different in the two coupling regions, changing from continuous (second-order) in the weak coupling limit to discontinuous (first-order) in the strong coupling limit. Estimates for the location of the multi-critical point, where the nature of the transition changes, have ranged from $`U_m1.5`$ to $`U_m5`$ (and $`V_mU_m/2`$). Despite the broad uncertainty in the actual value of the tricritical point, the phase diagram was believed to be well understood.

Recently, however, by studying the EHM ground state broken symmetries using level crossings in excitation spectra obtained by exact diagonalization, Nakamura  has argued for the existence of a novel bond-order-wave (BOW) phase for small to intermediate values of $`U`$ and $`V`$ in a narrow strip between the CDW and the SDW phases. The BOW phase is characterized by alternating strengths of the expectation value of the kinetic energy operator on the bonds. It is predicted to be a state where the discrete (two-fold) symmetry is broken and should hence exhibit true long-range order. Nakamura thus argues that the transition between CDW and SDW phases in this region is replaced two separate transitions: $`(i)`$ a continuous transition from CDW to BOW; and $`(ii)`$ a Kosterlitz-Thouless spin-gap transition from BOW to SDW. The BOW region vanishes at the multi-critical point beyond which the transition between CDW and SDW phases is direct and discontinuous. A schematic phase diagram including Nakamura’s BOW state is shown in Fig. 1.

Considering the long history of the 1D EHM and the large number of studies of the $`U2V`$ region with a variety of analytical and numerical tools, the proposal of a new phase is certainly remarkable. Importantly, the level crossing method used by Nakamura cannot by itself exclude the conventional scenario of a direct SDW-CDW transition for the whole range of $`U,V>0`$; a level crossing corresponding to this transition was also found  between the SDW-BOW and BOW-CDW crossing curves. The position of the BOW-CDW level crossing is, however, in closer agreement with the strong-coupling result for the vanishing of the CDW order, and this was taken as evidence of a long-range-ordered BOW in the ground state for certain parameters. It is important to confirm this hitherto undiscovered phase using other methods.

To attempt this confirmation, we have used the highly efficient stochastic series expansion (SSE) quantum Monte Carlo method  to study the EHM at half filling in the vicinity of $`U=2V`$. This method allows us to probe directly the spin- charge- and bond-order correlations in the ground state of lattices with more than one hundred sites (up to 256 sites were used in this study). Using finite-size scaling techniques for the various order parameters, we confirm the existence of a BOW state with spin and charge gaps in a region very close to that predicted by Nakamura for small $`U,V`$. We also further improved the SSE simulations by applying a quantum version of the thermal parallel tempering scheme (or exchange Monte Carlo)  for simulations close to and across the phase boundaries. This “quantum parallel tempering” greatly reduced the effects of “sticking” — where the simulation gets trapped in the wrong phase close to a phase boundary — and was found to be particularly useful for the discontinuous (first-order) direct SDW-CDW transition. As a consequence, we were able to obtain a more accurate estimate for the location of the multi-critical point $`(U_m,V_m)`$ where the BOW phase vanishes and is replaced by a first-order SDW-CDW transition line. As we discuss below, we find $`U_m=4.7\pm 0.1,V_m=2.51\pm 0.04`$.

In order to investigate the possibility of soliton lattices forming out of the long-range CDW and BOW states when doping away from half-filling, we have also carried out some simulations of lightly doped systems. We find that the in both cases the ground state is a Luther-Emery liquid, with a spin gap and slow algebraic decay ($`r^{K_\rho }`$, with $`K_\rho 0.5`$) of the dominant CDW and BOW correlations.

The remainder of the paper is organized into four sections and two appendices. In Sec. II we briefly sketch the SSE method and introduce the different observables we study. In Sec. III we present the results of our simulations at half-filling and discuss their interpretation. Doped systems are considered in Sec. IV. We conclude with a brief summary in Sec. V. In Appendix A we present some important details of the extension of the SSE method to allow efficient loop updates for fermions. We illustrate the advantages of the quantum parallel tempering scheme in Appendix B.

## II Numerical Methods and Observables

### A The SSE method and its fermion loop-update extension

The SSE method  is a finite-temperature quantum Monte Carlo method based on importance sampling of the diagonal elements of the Taylor expansion of $`e^{\beta H}`$, where $`\beta `$ is the inverse temperature $`\beta =t/T`$. Ground state expectation values can be obtained using sufficiently large values of $`\beta `$, and there are no approximations beyond statistical errors. Recently, in the context of spin systems, an efficient “operator loop update” was developed to sample the operator sequences appearing in the expansion. The resulting method has proven to be very efficient for several different models. To apply the most efficient variant of SSE method to the EHM, we need to generalize the previous operator loop update scheme to spinful fermions. This is an important extension, but because of its technical nature we have relegated our detailed discussion of it to an appendix.

We have applied the SSE algorithm to the 1D EHM for system sizes ranging from $`N=8`$ to $`256`$ sites, with maximum inverse temperatures $`\beta `$ chosen appropriately to isolate the ground state. We have verified the correctness of the simulation code by comparing $`N=8`$ results with exact diagonalization (Lanczos) results.

Although the operator-loop update is indeed significantly more efficient than previous local updates for sampling of the SSE configurations, we still have problems with “trapping” close to a first-order phase transition, i.e., the simulation can get stuck in the wrong phase very close to the critical point. There are also problems with slow dynamics in long-range ordered phases with a broken discrete symmetry (such as BOW or CDW phases). In order to overcome these problems we have developed a “quantum parallel tempering” scheme — a generalization of the thermal parallel tempering method  commonly used to equilibrate classical spin glass simulations. The method amounts to running several simulations on a parallel computer, using a fixed value of $`U`$ and different closely spaced values of $`V`$ at and around the critical value $`V_c`$. Along with the usual Monte Carlo updates, we attempt to swap the configurations for processes with adjacent values of $`V`$ at regular intervals (typically after every Monte Carlo step) according to a scheme that maintains detailed balance in the space of the parallel simulations, as explained in Appendix B. In contrast to Ref. , we here find parallel tempering to be particularly useful in the study of the first-order transition, where the problem of trapping is the most pronounced. In Appendix B we also present a comparative example to illustrate the improvement obtained by parallel tempering.

### B Observables

In addition to the ground state energy, $`E=H/N`$, the observables we study include the static structure factors and susceptibilities corresponding to the different phases (CDW, SDW and BOW). The structure factors are given by

$`S_{SDW}(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j,k}{}}e^{iq(jk)}S_j^zS_k^z,`$ (4)
$`S_{CDW}(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j,k}{}}e^{iq(jk)}n_jn_kn_j^2,`$ (5)
$`S_{BOW}(q)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{j,k}{}}e^{iq(jk)}k_jk_kk_j^2,`$ (6)
where

$$k_j=\underset{\sigma =,}{}(c_{j+1,\sigma }^{}c_{j,\sigma }+h.c.)$$
(7)
is the kinetic energy operator associated with the $`j`$th bond. The corresponding static susceptibilities are given by

$$\chi _{SDW}(q)=\frac{1}{N}\underset{j,k}{}e^{iq(jk)}_0^\beta 𝑑\tau S_j^z(\tau )S_k^z(0)$$
(8)
and analogous expressions for $`\chi _{CDW}(q)`$ and $`\chi _{BOW}(q)`$. Since all the phases mentioned have a period 2, the staggered structure factor and susceptibilities are the most important observables. We define order parameters for the phases in terms of the staggered structure factors:

$$m_\alpha =\sqrt{S_\alpha (\pi )/N},$$
(9)
where $`\alpha =`$ CDW, SDW, or BOW. We have also studied the charge stiffness constant, $`\rho _c`$. It is defined as the second derivative of the internal energy per site, $`E`$, with respect to a twist, $`\varphi `$,

$$\rho _c=\frac{^2E(\varphi )}{\varphi ^2},$$
(10)
under which the hopping term in the Hamiltonian (3) is replaced by

$$k_c(\varphi )=t\underset{j,\sigma }{}(e^{i\varphi }c_{j+1,\sigma }^{}c_{j,\sigma }+h.c.).$$
(11)
The spin stiffness constant, $`\rho _s`$, is defined by a similar expression, with the hopping term now being replaced by

$$k_s(\varphi )=t\underset{j,\sigma }{}(e^{i\varphi _\sigma }c_{j+1,\sigma }^{}c_{j,\sigma }+h.c.),$$
(12)
with $`\varphi _{}=\varphi _{}=\varphi `$. In the framework of the SSE method, the estimators for the charge and spin stiffness are given in terms of expectation values of squared winding numbers (see Appendix A).

## III Results at half-filling

As noted above, we have studied chains with $`N`$ up to $`256`$ with periodic boundary conditions. Typically, an inverse temperature of $`\beta =2N`$ was sufficient for the calculated properties to have converged to their ground state values, except in the case of $`N=256`$, for which $`\beta =4N`$ was needed for some quantities. In this section we first discuss our evidence for the existence of a long-range BOW phase, then our analysis of the continuous BOW-CDW and SDW-BOW transitions for small $`(U,V)`$, the discontinuous SDW-CDW transition for large $`(U,V)`$, and finally our determination of the location of the multi-critical point separating these transitions.

### A Existence of the BOW phase

Plots of the variation of the staggered susceptibilities corresponding to the three different phases — CDW, SDW, and BOW — show the existence of strong BOW fluctuations in a region with $`VU/2`$ in parameter space where Nakamura predicted a BOW state. Fig. 2 is such a plot for $`U=4`$ and $`1.7V<2.3`$. In a long-range ordered phase (BOW, CDW), the corresponding $`\chi (\pi )`$ is expected to diverge with increasing system whereas the other two susceptibilities should converge to constants. In the SDW phase there is no long-range order but algebraically decaying correlations of both SDW and BOW nature; hence $`\chi _{SDW}(\pi )`$ and $`\chi _{BOW}(\pi )`$ should both diverge here, but the BOW divergence should be much slower than in the long-ranged BOW phase. These behaviors are indeed seen in Fig. 2, with the susceptibilities for SDW, BOW, and CDW dominating in turn as $`V`$ is increased. The BOW-CDW phase boundary can be quite well resolved since it involves a standard second order (continuous) phase transition. On the other hand, the SDW-BOW boundary is more difficult to locate, for it involves a Kosterlitz-Thouless transition in which the spin gap opens exponentially slowly as one enters the BOW phase,  resulting in only a slow decay of the staggered SDW susceptibility in the BOW phase for the system sizes accessible in our work.

Fig. 3 shows ln$`[\chi _\alpha (\pi )]`$ and ln$`[S_\alpha (\pi )]`$ versus ln$`[N]`$ for the parameters $`(U,V)=(4,2.14)`$ for which the ground state should be inside the BOW phase. We find that both $`\chi _{BOW}(\pi )`$ and $`S_{BOW}(\pi )`$ diverge strongly with system size, whereas the structure factor and susceptibility corresponding to CDW have a maximum and then decrease with system size for large $`N`$. The SDW structure factor appears to have converged for $`N=256`$ but the susceptibility still shows a weak growth — in a spin-gapped BOW phase it should eventually converge, too, but if the gap is very small the convergence occurs only for much larger systems. The growth with $`N`$ seen here is much slower than $`N`$, which should be the asymptotic behavior in an SDW phase for any spin-rotationally invariant 1D system,  and the growth slows with increasing $`N`$. Hence an asymptotic divergence of $`\chi _{SDW}(\pi )`$ can be excluded. The dominant asymptotic characteristic of the ground state is clearly BOW. The system sizes considered are not large enough for $`S_{BOW}(\pi )`$ to have reached the asymptotic behavior $`N`$ expected if there is long-range order, which we will explain further below. The very fast divergence of $`\chi _{BOW}(\pi )`$ is expected on account of the two-fold degenerate BOW ground state. For finite $`N`$ this degeneracy is not perfect, but an exponentially fast closing of the gap between the symmetric and antisymmetric linear combinations of the two asymptotically degenerate symmetry-broken ordered states can be expected, which would eventually cause $`\chi _{BOW}(\pi )`$ to diverge exponentially.

The most direct evidence for a long-range BOW comes from the the real-space kinetic energy correlation function

$$C_{BOW}(r)=\frac{1}{N}\underset{i=1}{\overset{N}{}}k_ik_{i+r}k_i^2.$$
(13)
As seen in Figure 4, this correlation function oscillates with period 2 and its magnitude decays considerably for short distances. For long distances there is a convergence to a constant, non-zero magnitude, which is the same within statistical errors for $`N=128`$ and $`256`$. The significant enhancement of the correlations at short distances explain the deviations from the expected asymptotic linear scaling of the integrated correlation function, $`S_{BOW}(\pi )`$, for the system sizes shown in Figure 3.

Further proof of the existence of the BOW phase is obtained by looking for spin and charge gaps in this region. Instead of calculating the gaps directly, which can not easily be done to high accuracy for large system sizes, we use the following indirect method: It is known  that if the ground state of a 1D system is gapless in the spin sector, the Luttinger liquid parameter $`K_\sigma `$ governing the asymptotic equal-time spin correlation function is $`K_\sigma =1`$. It has been further shown  that the slope $`S_{SDW}(q)/q`$ gives $`K_\sigma /\pi `$ in the limit $`q0`$. Hence, $`S_{SDW}(q)/q1/\pi `$ as $`q0`$. On the other hand, if the ground state has a spin gap, $`S_{SDW}(q)/q0`$ as $`q0`$. With this criterion even a very small spin gap can be detected, since it is in practice sufficient to see that $`\pi S_{SDW}(q)/q`$ decays below $`1`$ for small $`q`$ to conclude that $`K_\sigma 1`$ and hence that a spin gap must be present. Similarly, for a ground state with no charge gap, $`\pi S_{CDW}(q)/qK_\rho `$ as $`q0`$ whereas if the ground state does have a charge gap, $`S_{CDW}(q)/q0`$ as $`q0`$. Unlike $`K_\sigma `$, where the value is fixed at $`1`$ for spin rotationally invariant systems, the Luttinger liquid charge correlation parameter $`K_\rho `$ is a function of $`U`$ and $`V`$, and its precise value for given $`U`$ and $`V`$ is not known (except at $`V=0`$ ). Due to the logarithmic corrections typical for 1D systems, it is very difficult to observe numerically that $`\pi S_{SDW}(q)/q`$ becomes exactly $`1`$. Empirically, we have found that in the gapless case the value $`1`$ is always approached from above (which is the case also for spin systems ), and hence the detection of the spin gap using this quantity is not hampered by the log-corrections — if $`\pi S_{SDW}(q)/q`$ decays below $`1`$ one can conclude that here is a gap.

Fig. 5 shows $`\pi S_{SDW}(q)/q`$ and $`\pi S_{CDW}(q)/q`$ versus $`q/\pi `$ for $`U=4`$ and two values of $`V`$. One of the points ($`V=2.14`$) is inside the BOW phase, whereas the other ($`V=1.8`$) is in the SDW phase. The $`\pi S_{SDW}(q)/q`$ curve for $`V=1.8`$ is close to $`1`$ for a wide range of $`q`$ values, whereas the $`V=2.14`$ curve exhibits a sharp drop as $`q0`$ indicating, respectively, the absence and the presence of a spin gap. Similarly, the evidence for a vanishing limit of $`S_{CDW}(q)/q`$ and hence of a charge gap for $`V=1.8`$ is clear. Since the point $`V=2.14`$ is quite close to the critical point ($`V_c=2.16`$) where the charge gap vanishes, the magnitude of the gap is very small and we need to go to still smaller $`q`$, i.e., larger system size, to see a pronounced effect like that for $`V=1.8`$. Nevertheless, the downturn for the smallest $`q`$ is a good indication of a gap.

The opening of spin and charge gaps can also be detected in the spin and charge stiffness constants, which should vanish as $`N\mathrm{}`$ if there are gaps. The asymptotic charge stiffness should hence be non-zero only exactly at the BOW-CDW phase boundary. The spin stiffness should be non-zero in the SDW phase, should approach a constant value exactly at the phase boundary (with logarithmic size-corrections), and vanish inside the CDW phase. In Fig. 6 we show the stiffness constants for $`U=4`$ in the neighborhood of the BOW phase. As expected, the charge stiffness peaks at the BOW-CDW phase boundary and decreases rapidly away from it, confirming the vanishing of the charge gap only at the phase boundary. The peak becomes very sharp for large system sizes, and the finite-size corrections to its location are small. We find this the most accurate way to locate the BOW-CDW phase boundary. The spin stiffness is clearly zero in the CDW phase, and a sharp decrease with increasing $`N`$ is also seen for $`V`$ values well inside the BOW phase. Since the spin gap opens up exponentially slowly at the SDW-BOW boundary it is difficult to locate the transition this way. Our data nevertheless indicate that the BOW phase at $`U=4`$ may not extend down to the value $`V1.82`$ obtained by Nakamura. We will discuss this phase transition and determine the transition point more accurately below, in Sec. III-C.

### B BOW-CDW transition

In addition to proving the existence of the BOW phase, we have studied in detail the nature of the continuous BOW-CDW transition for two different values of $`U`$ ($`U<U_m`$). For $`(U,V)=(U_c,V_c)`$, i.e., on the BOW-CDW phase boundary, the real space staggered charge and kinetic energy correlation functions fall off algebraically as

$`n_in_{i+r}(1)^r`$ $``$ $`r^\eta ,`$ (14)
$`(K_iK_{i+r}K_i^2)(1)^r`$ $``$ $`r^\eta .`$ (15)
Based on conformal field theory calculations for similar phase transitions in 1D spin systems, the exponent $`\eta `$ can be expected to depend on $`(U_c,V_c)`$ but should be the same for both the CDW and BOW correlations. This gives the finite-size scaling of the structure factor and the susceptibility at the critical point:

$`S_{CDW,BOW}(\pi )`$ $``$ $`N^{1\eta },`$ (16)
$`\chi _{CDW,BOW}(\pi )`$ $``$ $`N^{2\eta }.`$ (17)
With a spin gap but no charge gap, as was demonstrated above, we expect the critical state to be of the Luther-Emery liquid type. The exponent $`\eta `$ is then related to the Luttinger liquid parameter $`K_\rho `$ by $`\eta =1K_\rho `$.

Fig. 7 presents plots of ln$`[\chi _{CDW}]`$ and ln$`[\chi _{BOW}]`$ versus ln$`[N]`$ for $`U=4`$ and three different values of $`V`$ around the critical point, which as discussed above should be close to $`2.16`$. The data points for $`V`$=2.16 indeed fall almost on straight lines, indicating critical scaling for both the CDW and BOW fluctuations. The value of the critical exponent $`\eta `$, obtained from the slope of the $`V=2.16`$ curves for both $`\chi _{CDW}`$ and $`\chi _{BOW}`$ is $`\eta 0.5`$. The scaling of the structure factors, $`S_{CDW}`$ and $`S_{BOW}`$, at $`V=2.16`$ is also consistent with $`\eta 0.5`$. It is, however, difficult to extract a precise value for $`\eta `$ from this finite-size scaling, due to subleading corrections to the scaling, as well as effects from the fact that the $`U,V`$ point studied is not exactly on the phase boundary. As was discussed in Sec. III-A, the Luttinger liquid parameter $`K_\rho `$ can also be extracted from the $`q0`$ limit of $`S_{CDW}(q)/q`$. This is in general a more accurate method, since the convergence with system size is faster for the subleading $`1/r^2`$ contribution to the correlation function which this estimator accesses.  Fig. 8 shows results for $`U=4`$ and $`U=3`$ and the respective critical $`V`$-values. The $`q0`$ behavior gives $`K_\rho =0.44\pm 0.01`$ for $`U=4`$, i.e., $`\eta =0.56\pm 0.01`$, which hence is consistent with the finite-size scaling of the $`q=\pi `$ quantities. For $`U=3`$ we obtain $`V_c=1.65`$, in agreement with Nakamura’s result,  and the critical exponent $`\eta =0.47\pm 0.01`$.

### C SDW-BOW transition

The SDW-BOW transition is marked by the opening of a spin-gap in the electronic energy spectrum. As argued by Nakamura, it is a quantum phase transition of the Kosterlitz-Thouless type and therefore the gap opens up exponentially slowly. This makes it difficult to determine the phase boundary numerically. The numerical data is affected by large finite-size effects that persist up to very large system sizes. As discussed in Sec. IIIA, the most reliable evidence of the existence of a spin-gap is obtained from the behavior of $`S_{SDW}(q)/q`$ as $`q0`$. In practice, an asymptotic value of $`\pi S_{SDW}(q)/q<1`$ as $`q0`$ in any (large) system is an indication of the presence of a spin-gap in the thermodynamic limit. This allows us to detect the presence of very small spin gaps. Fig. 9 shows the behavior of $`\pi S_{SDW}(q)/q`$ for $`U=4`$ and different values of $`V`$. In the gapless region, logarithmic corrections  makes it difficult to observe the approach to $`1`$ as $`q0`$. In analogy with spin systems, we expect the leading log corrections to vanish at the point where the spin gap opens, and therefore exactly at the critical point there should be a clear scaling to $`1`$. An apparent reduction of the log correction is indeed seen in Fig. 9 as $`V`$ is increased towards $`1.88`$. Based on the results, we estimate the SDW-BOW boundary to be at $`V=1.89\pm 0.01`$ at $`U=4`$. This is slightly higher than Nakamura’s critical value $`V=1.82`$ for this $`U`$. We believe the difference is due to non-asymptotic finite-size effects in the exact diagonalization calculation, which used system sizes only up to $`N=14`$. Hence, we find that the BOW phase exists in a slightly smaller, while still significant, region of the phase space.

### D First-order SDW-CDW transition

For $`U>U_m`$, the transition is a discontinuous (first-order) direct SDW-CDW transition with no intervening BOW phase. Fig. 10 shows the $`V`$ dependence of the CDW order parameter, the total energy, and the kinetic energy across the phase boundary for $`U=8`$, which according to previous studies  should be well within the regime of first-order transitions. The characteristics of a first-order transition are indeed quite apparent. The order parameter and the kinetic energy change rapidly at the transition point, $`V_c4.14`$. The finite-size effects diminish with increasing $`N`$ as the results approach the limiting behavior of a discontinuity in the order parameter and the kinetic energy in the thermodynamic limit. The total energy remains continuous, but there is a clear break in slope at the transition.

The size dependence of the BOW order-parameter is shown in Fig. 11. It becomes considerably smaller inside the CDW phase than before the transition. This is expected, since in the SDW phase, but not in the CDW phase, there should be power-law decaying BOW correlations. The BOW order parameter decays rapidly with the system size, however, confirming that there is no long-range BOW for this $`U>U_m`$.

The behavior with increasingly sharp discontinuities seen in Fig. 10 and 11 indicates a first-order transition due to an avoided level crossing. Note that with increasing chain length the CDW order parameter approaches its thermodynamic value from above for $`V<V_c`$ and from below for $`V`$ just above $`V_c`$. The curves for different system sizes cross one another in the neighborhood of $`V=V_c`$ and then once again for a higher $`V`$. The second crossing point moves down towards the first one as $`N`$ increases, whereas the first crossing does not change much with $`V`$ and appears to be a good criterion for locating the transition point.

The two curve crossings can be understood as follows: In a transition caused by an avoided level crossing, a crossing of the order parameter curves close to the critical coupling (approaching the critical coupling as $`N\mathrm{}`$) can be expected since the low-energy levels corresponding to an ordered and disordered state swap characters within a a parameter range $`V\pm \mathrm{\Delta }_V(N)`$, with $`\mathrm{\Delta }_V(N)0`$ as $`N\mathrm{}`$. This behavior is seen clearly in Fig. 10. The finite-$`N`$ ground state starts to develop CDW characteristics at $`V\mathrm{\Delta }_V(N)`$ and thus, for a fixed $`V<V_c`$, the CDW order parameter decreases with increasing $`N`$. An analogous argument for fixed $`V>V_c`$ close to $`V_c`$ gives that the in this case the CDW order must increase with increasing $`N`$. On the other hand, for $`VV_c`$ the real-space CDW correlations are enhanced at short distances (in the same way as the BOW correlations shown in Fig. 4) and for small system sizes there is also some enhancement of the long-distance correlations due to the periodic boundary conditions. Hence, one can expect the CDW order parameter, when defined and measured in terms of its squared expectation Eq. (9), to again decrease with $`N`$ for $`VV_c`$ and this explains the second crossing of the order parameter curves seen in Fig. 10.

### E Multi-critical point

Although the existence of the tricritical point (which, in vew of the existence of the BOW phase, we refer to as the multi-critical point) separating the first-order and continuous transition to the CDW state has long been known, its location in the $`(U,V)`$ plane has not previously been determined accurately using large system sizes. Hirsch  estimated a value of $`U_m=3`$ using world line Monte Carlo. Cannon and Fradkin obtained $`U_m=1.5`$ using field theory techniques and world lines. Later Cannon, Scalettar and Fradkin obtained a value of $`U_m=3.55`$ using finite-size scaling of Lanczos results. Using a combination of bosonization and RG techniques, Voit  obtained $`U_m=4.76`$. However, as Voit also pointed out, the validity of bosonization and RG, which are applicable in the limit $`U,V0`$, for intermediate values of the parameters is a priori questionable.

By using larger system sizes and an alternative criterion to distinguish between a continuous transition and a first-order level crossing transition, we have obtained an estimate of the multi-critical point that we consider more accurate and reliable than the previous estimates. In contrast to most previous numerical studies, our method is not based on plotting histograms of the order parameter, although we will also present such histograms in the next section. In this section we first exploit the qualitatively different finite-size dependence of the growth of the order parameter close to the transition above and below the multi-critical point.

For fixed $`U`$, the order parameter curves for different system sizes cross each other at or very close to the critical point ($`V=V_c`$) in the case of a first-order transition, as discussed above in Sec. III-D. Such a crossing cannot occur at a continuous transition, where instead there should be finite-size scaling governed by Eq. (17). This qualitative difference in the finite-size dependence of the order parameter close to the transition point above and below the multi-critical point $`(U_m,V_m)`$ leads us to expect that in the neighborhood of this point, curves of the order parameter for different chain lengths will closely coincide with one another close to $`V=V_c`$, and $`U_m`$ is the point at which the curves barely touch each other. When the system size becomes sufficiently large one can also directly observe discontinuities developing when $`U>U_m`$, in the order parameter as well as in other quantities, as in Fig. 10. In practice this criterion, or any other criterion known to us, cannot be expected to be useful very close to the multi-critical point, where the transition is only weakly first-order and very large lattices are needed to detect discontinuities developing from avoided level crossings.

Fig. 12 shows the finite-size dependence of the CDW order parameter across the transition for three different values of $`U`$. For $`U=4.2`$, only the $`N=16`$ curve crosses the other curves, and this occurs far from the critical point (as determined using the peak in the charge stiffness, as discussed in Sec. III-A). The non-crossing for larger system sizes show that the transition must be continuous at this $`U`$. For $`U=5.2`$ all curves show crossing behavior and a discontinuity can also be seen developing for the largest system size, i.e., the transition is here of first order. The curves for $`U=4.6`$ closely follow the expected behavior at the multi-critical point, with the curves for the largest systems barely touching each other. Based on data also for other values of $`U`$ we estimate the multi-critical point to be $`(U_m=4.7\pm 0.1,V_m=2.51\pm 0.04)`$. This agrees very well with Voit’s estimate $`(U_m=4.76)`$. However, it is not clear whether this agreement is fortuitous or whether there is some underlying symmetry that renders bosonization and RG (that assume $`U,V1`$) applicable close to the multi-critical point.

### F CDW order parameter histograms

Previous studies of the multi-critical point have exploited the existence of a 3-peak structure in the distribution of the CDW order parameter for a discontinuous SDW-CDW transition in the vicinity of the critical point and its absence at a continuous transition. Outside the CDW phase, the distribution of the CDW order parameter is peaked around zero. For a continuous transition to a CDW state this peak splits into two (corresponding to positive and negative values of the order parameter), which gradually move apart from each other inside the CDW phase. In a first-order transition, on the other hand, the order parameter takes a non-zero value immediately as the CDW phase is entered and hence the two peaks emerge already separated from each other. Furthermore, at the phase boundary the CDW phase coexists with the competing phase, and this is reflected as a central peak remaining in the CDW order parameter distribution. The position of the multi-critical point can then in principle be obtained by locating the point where the 3-peak structure first appears. In practice, the accuracy of this method is limited by the fact that the discontinuity is very small for a first-order transition close to the multi-critical point and very large system sizes are then needed to observe the three peaks. This problem is, of course, common to all methods for distinguishing between a continuous and weakly first-order transition.

In his early QMC study, Hirsch observed a 3-peak structure even for $`U`$ as small as $`3`$ and therefore concluded that the transition is of first order already there. For larger $`U`$, an unexplained $`4`$-peak structure was seen. We have repeated histogram calculations for the lattice size $`N=32`$ studied by Hirsch. In Fig. 13 we show results for $`U=6`$, $`V=3.15`$, where a 4-peak structure was seen in the earlier calculation. We only find a central peak, which show that the system is not in the CDW state for these parameters. There are, however, already signs of side peaks developing, which shows that the system is close to the CDW phase. The significant differences with the earlier result could partially be errors due to the Trotter decomposition used in the world-line simulation method. Temperature effects are only minor, as also shown in Fig. 13. At $`\beta =8`$, which was used in Ref. , the histogram is only slightly more sharply peaked than at $`\beta =16`$ and $`32`$. Most likely, the simulation giving the $`4`$-peak structure was not sufficiently long, as it consisted of only $`10^4`$ Monte Carlo steps. Even with the more efficient SSE algorithm used in the present work, we find that the autocorrelation times are quite long close to the first order transition (see Appendix B) and short simulation can produce incorrect order parameter histograms similar to those shown in Ref. . For the histograms shown here, on the order of $`10^710^8`$ SSE Monte Carlo steps were used.

In Fig. 13 we also show results for several values of $`V`$ across the phase transition. A clear 3-peak structure (i.e., three peaks in the range $`m_{CDW}[1,1]`$, of which we only show the positive part) with peaks of almost the same heights can be seen for $`V=3.165`$. In Fig. 14 we show results for $`N=64`$. At $`U=6`$, the 3-peak structure appears for $`V3.156`$, i.e, at a value slightly lower than for the $`N=32`$ system. The size of the $`V`$ region in which three peaks can be observed is also significantly smaller, reflecting the sharpening of the first order transition cause by an avoided level crossing. At $`U=5`$, which we have argued above should be close to but above the multi-critical point, we do not observe three peaks. The histogram becomes very flat for an extended range of $`m_{CDW}`$, however, and the side peak emerges at a finite value of $`m_{CDW}`$. This is consistent with the transition still being of first order at $`V=5`$. Going to still lower $`V`$-values, the peak just becomes narrower, and it is not possible to definitely conclude this way when the transition becomes continuous.

## IV Doped systems

An interesting question naturally arises from the existence of the Nakamura BOW phase: Can the EHM model support a soliton lattice when doped slightly away from half-filling? Such a state exists in the Su-Schrieffer-Heeger model in the adiabatic limit (with classical, “frozen” phonons),  and also when electron-electron interactions are taken into account. In these models, the quantum nature of the phonon field is not taken into account, however. It is known that the dimerized state at half-filling survives even in the presence of quantum fluctuations, at least up to a critical value of the phonon frequency. However, to our knowledge, there have been no reliable numerical calculations addressing the stability of the soliton lattice in the presence of fully quantum mechanical phonons. The Nakamura BOW is similar to a dimerized lattice with quantum fluctuations, and hence a study of its evolution with hole doping can give insights also into the quantum phonon problem. There are also unresolved issues regarding the doped CDW state. In order to investigate the evolution of the long-range ordered states upon doping, we have studied the EHM model also away from half-filling, focusing on two parameter values in which the half-filled system is in the BOW (using $`U=4`$, $`V=2.14`$) or CDW phase (using $`U=4`$, $`V=2.5`$).

We first discuss the effects of doping on the spin and charge gaps in the half-filled CDW and BOW ground states. As in previous sections, we make use of the behavior of the static structure factor in the limit $`q0`$. Figs. 15 and 16 show $`\pi S_{CDW}(q)/q`$ and $`\pi S_{SDW}(q)/q`$ as a function of $`q`$ for a range of doping levels, both for parameters where the half-filled system is in the BOW phase (Fig. 15) and in the CDW phase (Fig. 16). From the data we conclude that upon doping away from half-filling, the charge gap vanishes immediately whereas the spin gap survives. This is true for both the CDW and BOW parent states. This behavior is characteristic of a Luther-Emery liquid,  in which the charge sector can be described in terms of a Luttinger liquid and the spin sector is gapped. The limiting value of $`\pi S_c(q)/q`$ as $`q0`$ indicates that the Luttinger liquid exponent $`K_\rho 0.5`$ in both the cases, with only a weak dependence on the doping level for the parameters considered here. Note the cross-over behavior occurring in the charge structure at $`q2\pi \delta =4k_F`$ in Fig. 16 (which is accompanied by a peak in the corresponding susceptibility, as will be shown below), reflecting a weak repulsion between dopant holes. No cross-over in the charge structure is seen in Fig. 15, where the parent state is a BOW.

Fig. 17 shows the variation of the ground state static susceptibilities for several doping levels in a chain of length $`N=128`$ for the parameters $`U=4,V=2.5`$. For $`\delta >0`$ the charge susceptibility converges to a non-zero value as $`q0`$, again showing the absence of a charge gap. Very strong $`2k_F`$ peaks are evident, and weaker $`4k_F`$ peaks are also clearly visible. The $`2k_F`$ peaks diverge with the system size whereas the $`4k_F`$ peaks are non-divergent, in accord with the Luther-Emery picture. For a Luther-Emery liquid, the charge correlations decay with distance $`r`$ as $`r^{K_\rho }`$, which gives $`\chi _{CDW}(2k_F)N^{2K_\rho }`$ for the finite-size scaling of the corresponding $`2k_F`$ susceptibility. Fig. 18 shows the size dependence for $`\delta =0.0625`$ on a log-log scale. For system sizes $`N64`$, the data is seen to fall on a line with slope $`1.5`$, consistent with the value $`K_\rho 0.5`$ extracted above.

As a further test of the Luther-Emery liquid nature of the ground state away from half-filling, we have studied the real-space charge and bond correlations as a function of distance. Fig. 19 shows the charge correlation for two different system sizes at a dopant concentration of 6.25% and interaction parameters $`U=4`$ and $`V=2.5`$. The ground state at half-filling is a CDW, and for the doped system we find solitonic features with alternating A and B phases separated by domain walls or kinks. The correlation decay with distance, however, and there is no real soliton lattice. In fact, the decay of the magnitude of the peaks is well approximated by an envelope curve of the form $`yx^{0.5}`$, and hence also these data are consistent with a Luther-Emery state with $`K_\rho 0.5`$.

Fig. 20 shows a similar plot of the real-space bond-order correlation for $`U=4`$ and $`V=2.14`$ at the same dopant concentration and for the same system sizes. The ground state of the half-filled system for this choice of parameters is here a BOW, and away from half-filling, the dominant correlation are still of bond-order type. Once again, the ground state of the doped system has solitonic features with an algebraic decay of the magnitude of the peaks. As in the previous case, the decay is consistent with a Luther-Emery liquid with $`K_\rho 0.5`$.

## V Summary

To summarize, we have studied the 1D EHM using the SSE method incorporating an efficient operator-loop update and a “quantum parallel tempering” scheme. Our results confirm the surprising prediction  of the existence of a novel long-range-ordered BOW phase between the well-known CDW and SDW phases in the ground state phase diagram for small to intermediate values of the on-site interaction $`U`$ ($`U<U_m`$). We have presented several ways to detect the spin and charge gaps expected in the BOW phase and have also probed directly the BOW correlations and concluded that true long-range order develops. We have studied a few points on the BOW-CDW phase boundary and obtained a very good agreement with Nakamura’s level crossing prediction  for the location of this phase boundary. For the SDW-BOW phase boundary, our results indicate a higher critical $`V`$ for fixed $`U`$ than given by the level crossing method and thus over-all a slightly smaller size of the BOW phase. Our results are for significantly larger systems than in the previous study and it is not surprising that the finite-size effects in the level crossings can be large for the SDW-BOW transition since the spin gap opens exponentially slowly in this Kosterlitz-Thouless transition. An over-estimation of the size of the BOW phase from the level crossings is also apparent considering that our estimated multi-critical point is well within the BOW phase of Nakamura’s phase digram. Since our BOW-CDW phase boundaries agree, this indicates problems with the scaling of the exact SDW-BOW level crossings close to the multi-critical point, as was also mentioned by Nakamura. For large values of $`U`$ ($`U>U_m`$) the transition is discontinuous (first-order). We have shown that curves of the CDW order parameter across this boundary for different system sizes cross each other twice, and explained this behavior in terms of an avoided level crossing. We have also used the curve crossings as a means to locate the position of the multi-critical point with greater accuracy than previously attained. Our estimate for the multi-critical point is $`U_m=4.7\pm 0.1,V_m=2.51\pm 0.04`$.

We have also studied systems doped slightly away from half-filling. We find that both the doped CDW and BOW states give rise to ground states of the Luther-Emery type, i.e., the quantum fluctuations do not allow the formation of true soliton lattices. Based on the fact that the BOW state is very similar to the dimerized ground state of models with finite-frequency (non-adiabatic) phonons, we conjecture that the soliton lattice is also unstable to arbitrarily weak quantum fluctuations in these models, unless two- or three-dimensional couplings are taken into account.

After our completion of the numerical calculations at half-filling Tsuchiizu and Furusaki  presented a weak-coupling g-ology calculation taking into account second-order corrections to the coupling constants. They obtained a phase diagram in very good quantitative agreement with ours, including the location of the multi-critical point.

###### Acknowledgements.

We would like to thank R. T. Clay for discussions of the operator-loop update for fermions. We thank M. Nakamura for sending us some of his numerical results for comparisons. A.W.S. would like to thank O. Sushkov for discussions. This work was supported by the NSF under grant No. DMR-97-12765 and by the Academy of Finland (project 26175). Most of the numerical calculations were carried out on the SGI Origin2000 and Condor systems at the NCSA, Urbana, Illinois. Some simulations were also carried out on the Origin2000 at CSC - Scientific Computing Ltd. in Helsinki.

## A Operator-loop updates in the SSE method.

The basic SSE approach has been discussed in several papers. . We here start with a brief review as a basis for introducing the operator-loop update  in the context of fermion models.

To implement the SSE method, the Hamiltonian (3) is written, upto an additive constant, in the form

$$H=\underset{b=1}{\overset{N}{}}(H_{1,b}+H_{2,b}+H_{3,b}),$$
(A1)
where $`b`$ is the bond connecting the sites $`b`$ and $`b+1`$, $`N`$ is the length of the chain, and the operators $`H_{a,b},a=1,2,3`$ are defined as

$`H_{1,b}`$ $`=`$ $`C\frac{U}{2}(n_{b,}\frac{1}{2})(n_{b,}\frac{1}{2})`$ (A2)
$``$ $`\frac{U}{2}(n_{b+1,}\frac{1}{2})(n_{b+1,}\frac{1}{2})`$ (A3)
$``$ $`V(n_b1)(n_{b+1}1),`$ (A4)
$`H_{2,b}`$ $`=`$ $`t(c_{b+1,}^{}c_{b,}+h.c.),`$ (A5)
$`H_{3,b}`$ $`=`$ $`t(c_{b+1,}^{}c_{b,}+h.c.).`$ (A6)
The constant $`C`$ shifts the zero of the energy and is chosen to ensure a non-negative expectation value for $`H_{1,b}`$ (needed in order to ensure a positive definite expansion of the partition function). Introducing a basis $`\{|\alpha \}=\{|\zeta _1,\zeta _2,\mathrm{},\zeta _N\}`$, where $`\zeta _i\{0,,,\}`$ denotes the electron state at the site $`i`$, the partition function $`Z`$=Tr$`\{e^{\beta H}\}`$ can be expanded in a Taylor series as

$$Z=\underset{\alpha }{}\underset{n=0}{\overset{\mathrm{}}{}}\underset{S_n}{}\frac{\beta ^n}{n!}\alpha \left|\underset{p=1}{\overset{n}{}}H_{a_p,b_p}\right|\alpha ,$$
(A7)
where $`S_n`$ denotes a sequence of index pairs defining the operator string $`_{p=1}^nH_{a_p,b_p}`$:

$$S_n=[a,b]_1[a,b]_2\mathrm{}[a,b]_n,$$
(A8)
where we use the notation $`[a,b]_p=[a_p,b_p]`$ and $`a\{1,2,3\}`$, $`b\{1,\mathrm{},N\}`$. In order to construct an efficient updating scheme, the Taylor series is truncated at a self-consistently determined power $`L`$, large enough to cause only an exponentially small, completely negligible error ($`L\beta |E|`$, where $`E`$ is the total internal energy; for details see Refs. ). We can then define a sampling space where the length of the sequences is fixed, by inserting $`Ln`$ unit operators, denoted by $`H_{0,0}`$, into each sequence. The terms in the partition function must be divided by $`\left(\genfrac{}{}{0pt}{}{L}{n}\right)`$ in order to compensate for the different ways of inserting the unit operators. The summation over $`n`$ is then implicitly included in the summation over all sequences of length $`L`$. The partition function takes the form

$$Z=\underset{\alpha }{}\underset{S_L}{}\frac{\beta ^n(Ln)!}{L!}\alpha \left|\underset{p=1}{\overset{L}{}}H_{a_p,b_p}\right|\alpha ,$$
(A9)
where the operator-index pairs $`[a,b]_p`$ now have $`a\{1,2,3\}`$ and $`b\{1,\mathrm{},N\}`$ or $`[a,b]_p=[0,0]`$. For convenience, we introduce a notation for states obtained by the action of the first $`p`$ elements of the operator string $`S_L`$:

$$|\alpha (p)\underset{j=1}{\overset{p}{}}H_{a_j,b_j}|\alpha .$$
(A10)
For a nonzero contribution to the partition function, $`|\alpha (L)=|\alpha (0)`$.

A Monte Carlo scheme is used to sample the configurations $`(\alpha ,S_L)`$ according to their relative contributions (weights) to $`Z`$. The sampling scheme consists of two types of updates, referred to as diagonal update and operator-loop updates. The diagonal update involves local substitutions of the form $`[0,0]_p[1,b]_p`$ and is attempted consecutively for every $`p\{1,\mathrm{},L\}`$ in the sequence for which $`[a,b]_p=[0,0]_p`$ or $`[1,b]_p`$. The updates are accepted with probabilities

$`P([0,0]_p[1,b]_p)`$ $`=`$ $`{\displaystyle \frac{N\beta M_{1,b}(p)}{Ln}},`$ (A11)
$`P([1,b]_p[0,0]_p)`$ $`=`$ $`{\displaystyle \frac{Ln+1}{N\beta M_{1,b}(p)}},`$ (A12)
where

$$M_{a,b}(p)=\zeta _b(p),\zeta _{b+1}(p)|H_{a,b}|\zeta _b(p1),\zeta _{b+1}(p1)$$
(A13)
is a matrix element on bond $`b`$, which in this case is diagonal ($`a=1`$). Only a single state $`|\alpha (p)`$ is stored in the computer during the diagonal update. When off-diagonal operators are encountered during the successive scanning of the operator string, the corresponding electron states are updated so that the information needed for evaluation of the probabilities (A12) is always available when needed.

The operator-loop update has been discussed in detail in Ref.  in the context of spins. Here we present the construction of loops for fermions. As explained in Ref. , the matrix element in Eq. (A9) can be graphically represented by a set of $`n`$ vertices (corresponding to the $`n`$ non-unit operators in $`S_L`$) connected to one another by the propagated electron states. Each vertex has four “legs” with electron states $`|\zeta _i(p1),\zeta _{i+1}(p1)`$ and $`|\zeta _i(p),\zeta _{i+1}(p)`$ before and after the action of the associated Hamiltonian operator $`H_{a_p,b_p}`$. There are 32 allowed vertices – 16 diagonal ones and 8 each associated with the off-diagonal $`H_{2,b}`$ and $`H_{3,b}`$ \[see Fig. 21(a)\]. A configuration $`(\alpha ,S_L)`$ is completely specified by the leg states of the $`n`$ vertices — except for sites that do not have any operators acting on them.

To carry out the operator-loop update, the linked list of the $`n`$ vertices is first constructed. In addition to the electron states at the legs of each vertex, the list also contains the addresses (i.e., the location in $`S_L`$) of the next vertex and the corresponding leg that each leg is connected to. The loop construction begins with randomly choosing a vertex and an “entry” leg. The electron state at the entry leg is changed to one of the 3 other allowed states chosen at random. Next an “exit” leg is chosen (following a procedure described below) and its associated electron state is updated so that the new leg states constitute an allowed vertex \[see Fig. 21(b)\]. The exit leg will be linked to a leg of another vertex (or, if there is only one operator in the configuration which acts on the site in question, another leg on the same vertex) and this will be the entry leg for the next vertex. The electron state at this new entry leg is then updated to match the state at the exit leg of the previous vertex. A new exit leg is then chosen following the same procedure. This is repeated until the exit leg from a vertex points to the starting point of the loop, which implies that the loop is closed and a new allowed configuration has been generated.

To choose an exit leg — given a vertex, an entry leg and the updated electron state at the entry leg — all the legs can be considered in turn and attempts made to update the associated electron state so that the new leg states constitute an allowed vertex. Because of spin an charge conservation on the vertices, at a given exit leg there is at most one possible update of the electron state that can lead to an allowed vertex. Hence, the exit leg uniquely determines the new vertex and the probability of choosing a given leg should be proportional to the weight of the new vertex, i.e., a matrix element of the form (A13), which in this case can be either diagonal or off-diagonal. In practice, a fast selection of an exit leg and updating of the vertex state is achieved using two pre-generated tables. The first one contains the cumulative probabilities of the 4 exit legs given an entrance leg, the old vertex state, and the new state at the entrance. The second table contains the new vertex states corresponding to the updated entrance and exit legs.

A special case occurs if the initial update at the entry leg of the first vertex of a loop is a spin-flip, i.e. the electron state changes from $``$ to $``$ or vice versa. In this case, the vertex weight does not change when updated and as a consequence the “bounce process”, where the exit leg is the same as the entrance leg, does not have to be included in the loop construction. The loop then becomes deterministic, i.e., there is a unique exit leg given the entrance leg. This is similar to the “loop-exchange” algorithm proposed in the context of the world-line method.

A full Monte Carlo updating cycle (MC step) consists of a diagonal update, followed by the construction of a linked vertex list. Next a number of operator-loop updates are carried out and finally the vertices are mapped back into a corresponding sequence $`S_L`$. The loop update typically also implies changes in the stored state $`|\alpha =|\alpha (0)`$, as some of the vertex legs (links) span across the periodic boundary in the propagation direction. The number of up and down electrons can be changed by the operator-loop update, as can the spatial winding numbers, and the algorithm is hence fully grand canonical. Note that at high and moderately low temperatures there are typically some sites of the system which has no vertices associated with them. The states on these sites can be randomly changed, since they have no affect the configuration weight.

The number of loops constructed for every MC step is determined such that on an average a total of $`L`$ vertices (we typically use $`2L`$) are visited. The truncation $`L`$ and the number of loops are adjusted during he equilibration part of the simulation and are thereafter held fixed. $`L`$ is determined by requiring that the highest $`n`$ reached during equilibration is at most $`7080\%`$ of $`L`$.

In certain parameter regions the length of a loop can sometimes become extremely long before it closes — in practice, it may even never close. It is therefore necessary to impose a maximum length, beyond which the loop construction is terminated and a new starting point is chosen (typically, we use $`50L`$ for this cut-off length). In order to reduce the likelihood of the next loop also exceeding the termination length it has proven useful to carry out a diagonal update before starting the next loop. The loop termination does not violate detailed balance and does not cause any systematical errors in the results. In most cases, incomplete loop termination occurs so infrequently that it does not adversely affect the simulation. In analogy with Ref. , where a scheme (there called “worm” update) similar to the operator-loops considered here was first introduced within the continuous worl-line representation, the end points of the loop during construction can be related to the single-particle Green’s function of the system and hence the tendency for loops to become exceedingly long for some parameter values must be related to some physical properties of the system. This issue should be studied further.

Estimators for the various structure factors and susceptibilities have been discussed in previous articles. Here we only note that the charge and spin stiffness constants, Eq. (10), can be expressed in terms of spin and charge current operators in analogy with the spin stiffness of the Heisenberg antiferromagnet previously discussed in Ref. , leading to

$$\rho _{c,s}=\frac{[(n_R^{}n_L^{})\pm (n_R^{}n_L^{})]^2}{N\beta },$$
(A14)
where $`n_{R,L}^\sigma `$ are the number of kinetic energy operators in the SSE term propagating spin-$`\sigma `$ particles in the “right” and “left” direction on the ring. Because of spin and charge conservation, the topological winding numbers $`(n_R^\sigma n_L^\sigma )/N`$ can take only integer values.

Although the operator-loop algorithm very significantly speeds up SSE simulations, in many cases reducing the autocorrelation function by orders of magnitude, the dynamics is still very slow in some parameter regions. For the extended Hubbard model studied here, problems with very long autocorrelation times occur in the long-range ordered BOW and CDW phases. The problems are particularly severe for large systems close to the BOW-CDW phase boundary, where “trapping”in the wrong phase often occurs. The slow dynamics in the BOW phase is illustrated in Figure 22, which shows the simulation time dependence of $`S_{BOW}(\pi )`$ and $`\chi _{BOW}(\pi )`$ during a simulation of a $`256`$-site system at $`\beta =512`$ \[$`S_{BOW}`$ has converged at this $`\beta `$ but $`\chi _{BOW}`$ is about $`20\%`$ larger still at $`\beta =1024`$\]. It is evident that the BOW autocorrelation time here is tens of thousands of MC steps. The BOW susceptibility exhibits a behavior where it sometimes takes very small values (less that $`10^3`$ of the average value), but corresponding large fluctuations upwards do not occur, i.e., the distribution of the $`\chi _{BOW}(\pi )`$ estimator for individual configurations is very skewed. The structure factor exhibits a more symmetric distribution. This behavior can be understood as a consequence of the BOW ground state for a finite system being a symmetric combination of the two possible real-space symmetry-broken states. The symmetry is not broken in a finite system and the simulation is also not trapped in one of the real space-states. Hence, the wave function that is sampled in the simulation contains both the real-space states and the behavior seen in Fig. 22 indicates that individual configurations also contain both components, in such a way that transitions (“tunneling”) between the two real-space states can occur during the SSE propagation (which can be simply related  to a propagation in imaginary time), at least for some configurations. Tunneling can be inferred from the qualitatively different evolutions of the structure factor and the susceptibility in MC time. The susceptibility is an integral of the bond-order correlation, as in Eq. (8), which in configurations where tunneling occurs can be much smaller than in configurations with no tunneling, since correlations between states with the same real-space configuration contributes positively but correlations between different states give a negative contribution. The structure factor, on the other hand, is an equal time correlation function and would not be much reduced by tunneling if the tunneling times are short. This explains the qualitatively different distributions of the $`\chi _{BOW}`$ and $`S_{BOW}`$ measurements in Fig. 22. Evidently, the updating process is very slow in adding and removing tunneling events in the configurations, which maybe is not that surprising considering that the tunneling is between two states with a discrete broken symmetry. These problems do not occur in SSE simulations of systems with a broken continuous symmetry, such as the two-dimensional Heisenberg model.

The trapping and tunneling problems can be significantly reduced by using the parallel tempering scheme (or exchange Monte Carlo),  which is discussed below in Appendix B.

## B Quantum parallel tempering

The “quantum parallel tempering” scheme is a straight-forward generalization of the thermal parallel tempering  method commonly used to equilibrate classical spin glass simulations. Our implementation amounts to running several simulations simultaneously on a parallel computer, using a a fixed value of $`U`$ and different closely spaced values of $`V`$. Along with the usual Monte Carlo updates, we attempt to swap the configurations for processes with adjacent values of $`V`$ at regular intervals, typically after every Monte Carlo step, according to a scheme that maintains detailed balance in the extended ensemble of parallel simulations. The probability of swapping the $`V`$-values of runs $`i`$ and $`i+1`$, which are running at $`V_i`$ and $`V_{i+1}`$, respectively, before the swap, is

$$P_{\mathrm{swap}}(V_i,V_{i+1})=\mathrm{min}[1,\frac{\mathrm{W}_\mathrm{i}(\mathrm{V}_{\mathrm{i}+1})\mathrm{W}_{\mathrm{i}+1}(\mathrm{V}_\mathrm{i})}{\mathrm{W}_\mathrm{i}(\mathrm{V}_\mathrm{i})\mathrm{W}_{\mathrm{i}+1}(\mathrm{V}_{\mathrm{i}+1})}],$$
(B1)
where $`W_i(V)`$ is the weight of the $`i`$th simulation configuration evaluated with the coupling $`V`$. The swap probabilities for fixed $`\mathrm{\Delta }V=V_{i+1}V_i`$ decreases with increasing system size and decreasing temperature and hence $`\mathrm{\Delta }V`$ and the range of $`V`$-values (if the number of processes is fixed) must be chosen smaller for larger system sizes.

The computational effort required for the swapping process is very minor compared to the actual quantum Monte Carlo simulations. It is therefore useful to carry out several swap attempts of all pairs of neighboring simulations between every MC step. Histograms containing the number of times each of the current configurations has “occupied” each $`V`$-bin can then be constructed and used for adding the contributions of each configuration to all the $`V`$-bins. This can contribute to reducing the statistical error of measured quantities.

To illustrate the advantage of quantum parallel tempering, we show two sets of data — obtained with and without the use of tempering — for a system undergoing a first-order transition. Fig. 23 shows the CDW order parameter across the first-order SDW-CDW phase boundary at $`U=8`$. The upper panel shows the data obtained from individual runs; the lower panel shows data for the same parameters obtained using tempering. The length of the individual simulations was $`10^5`$ MC steps for all $`V`$ values, and this was also the number of steps performed by each process in the tempering runs. The improvement in the quality of the tempering data is evident, especially close to the transition point where two of the individual simulations have relaxed into the wrong phases. The statistical errors are hence severely underestimated due to the failure to equilibrate properly within the simulation time. The tempering error bars are also large at the phase transition, but in contrast to those of the individual simulations they are accurate error estimates. The errors rapidly become much smaller as one moves away from the transition point. The effects of tempering are also favorable further inside the CDW phase, where several of the individual simulations are apparently affected by trapping in configurations with defects, where the order is reduced.

The tempering acceptance rate during the run spanning across the phase transition in Fig. 23 is shown in Fig. 24. There is a sharp reduction in the acceptance rate at the transition. This reflects the rapid change in the SSE configurations across the phase boundary, which implies that the configuration weights evaluated with $`V`$ values from the “wrong” phase are likely to decrease and the swap according to the probability (B1) to be rejected.

Finally, we note that tempering, in general, is an application where a superlinear speed-up can be achieved in practice on parallel computers. In addition to doubling the density of data points when the number of processes is doubled, the statistical errors are also reduced. Sometimes the error reduction can be dramatic, but even in cases where there are no real problems with the dynamics of individual simulations the effects of tempering are often very favorable.


format: text

text:
# A condition on the chiral symmetry breaking solution of the Dyson-Schwinger equation in three-dimensional QED

(13,Apr 2001)

## Abstract

In three-dimensional QED, which is analyzed in the 1/$`N`$ expansion, we obtain a sufficient and necessary condition for a nontrivial solution of the Dyson-Schwinger equation to be chiral symmetry breaking solution. In the derivation, a normalization condition of the Goldstone bound state is used. It is showed that the existent analytical solutions satisfy this condition.

The interest in quantum electrodynamics in two space and one time dimensions (QED3) is twofold. On the one hand, QED3 itself as a quantum field theory exhibits many interesting features, such as dynamical chiral symmetry breaking (DCSB) \[1-5\] and confinement , which also exist in QCD4 and other gauge theories. Furthermore, unlike four dimensional QED, QED3 has a dimensional coupling constant $`e`$ and therefore is a superrenormalizable field theory, which makes it easier to be treated. It is expected that the study on QED3 would shed light on our understandings of more complicated gauge theories like QCD. On the other hand, QED3 has been used to describe some planar condensed matter systems, especially high temperature cuprates. It is well established that an undoped cuprate is a 2d quantum antiferromagnet and described by the spin-$`1/2`$ Heisenberg model , which can be mapped into a theory of massless fermions coupled to a U(1) gauge field . The chiral symmetry breaking in this U(1) gauge theory corresponds to the Néel ordering in the antiferromagnet .

One of the most interesting characteristics of QED3 is that it may exhibit DCSB, from which a massless fermion acquires a dynamically generated mass. This is a nonperturbative effect and conventional perturbative expansion approach is unable to investigate this problem. The standard mathematical tool of analyzing chiral phase transition is the Dyson-Schwinger (DS) nonlinear integral equation for the fermion self-energy. If the DS equation has only vanishing solution, the theory is chirally symmetric and the fermions remain massless. The case of interest to us is when a nontrivial solution of the DS equation occurs. Generally, the existence of a nontrivial solution of the DS equation was believed to lead unambiguously to the existence of DCSB. However, as pointed out by Cheng and Kuo , where they discussed this problem in quenched planar QED4, this is not always the case. Actually, the breaking of a global chiral symmetry is always accompanied by a (only one in U(1) gauge theory) massless Goldstone boson , which is a pseudoscalar bound state composed of a fermion and an antifermion . Once the DS equation develops a nontrivial solution, there should also be a nontrivial solution for the corresponding Bethe-Salpeter (BS) equation satisfied by this bound state wave function. Note that: these solutions are not independent and can be converted into each other under a simple transformation, and the BS wave function must satisfy an additional normalization condition to ensure the bound state is stable . Using these two facts, we may obtain a constraint on the nontrivial solution of the DS equation from the normalization condition of the BS wave function. Therefore, that the DS equation has a nontrivial solution is only a necessary but not a sufficient condition to lead to DCSB.

In the spirit of the above discussions, Cheng and Kuo  have considered the quenched planar QED4 and acquired a necessary and sufficient condition for a nontrivial solution of the DS equation to be symmetry breaking. They applied this condition to the explicit solution of the DS equation given by Kondo $`et`$ $`al.`$ , and found that, although DCSB occurs when there is a cutoff, once the cutoff is taken away DCSB disappears.

In this paper, we examine this problem in the continuous QED3 and verify whether or not the existed nontrivial solutions of the DS equations lead to DCSB. The Lagrangian of QED3 in Euclidean space is

$$L=\underset{i=1}{\overset{N}{}}\overline{\psi }_i(i_\mu +eA_\mu )\gamma _\mu \psi _i+\frac{1}{4}F_{\mu \nu }^2,$$
(1)
with $`N`$ fermions which are four-component spinors. In this case, the 4$`\times 4`$ $`\gamma _3`$ and $`\gamma _5`$ matrices can be well constructed, which anticommute with $`\gamma _0,`$ $`\gamma _1`$, and $`\gamma _2`$. For details about the $`\gamma `$ matrices, see Ref. . The Lagrangian is invariant under the chiral transformations $`\psi \mathrm{exp}(i\theta \gamma _3)\psi `$ and $`\psi \mathrm{exp}(i\omega \gamma _5)\psi `$ because the fermions are massless. A mass term $`m\overline{\psi }\psi ,`$ no matter it is added to the Lagrangian by hand or dynamically generated, would break these symmetries. Since the mechanism of dynamical mass generation for fermions needs no additional Higgs particles, it is very interesting to study DCSB in various field theories.

Power-counting arguments show that the square of the coupling constant $`e^2`$ has dimensions of mass, therefore QED3 is a superrenormalizable field theory and ultraviolet (UV) finite. However, while the UV behavior becomes better, in the massless case, perturbative expansions in $`e^2`$ lead to infrared (IR) divergences which appear already at two loops. One way to overcome this difficulty is using the 1/$`N`$ expansion .

We write the full fermion propagator as

$$S^1(p)=i\gamma pA(p^2)+B(p^2),$$
(2)
where $`A(p^2)`$ is the wave-function renormalization and $`B(p^2)`$ is the fermion self-energy. The DS equation for the full fermion propagator  is

$$S^1(p)=S_0^1(p)e^2\frac{d^3q}{(2\pi )^3}\gamma _\mu S(q)\mathrm{\Gamma }_\nu (p,q)D_{\mu \nu }(pq),$$
(3)
where $`\mathrm{\Gamma }_\nu (p,q)`$ is the full vertex function and the full photon propagator is

$$D_{\mu \nu }(q)=\left(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2}\right)\frac{1}{q^2[1+\mathrm{\Pi }(q^2)]}+\xi \frac{q_\mu q_\nu }{q^4},$$
(4)
with the vacuum polarization $`\mathrm{\Pi }(q^2)`$ related to the vacuum polarization tensor $`\mathrm{\Pi }_{\mu \nu }(q)`$ as follows:

$$\mathrm{\Pi }_{\mu \nu }(q)=(q^2\delta _{\mu \nu }q_\mu q_\nu )\mathrm{\Pi }(q^2).$$
(5)
Substituting Eq.(2) into Eq.(3) and taking trace on both sides, we obtain the equation for $`B(p^2)`$

$$B(p^2)=\frac{e^2}{4}\frac{d^3q}{(2\pi )^3}tr[\gamma _\mu S(q)\mathrm{\Gamma }_\nu (p,q)D_{\mu \nu }(pq)].$$
(6)
Multiplying both sides of Eq.(3) by $`\gamma p`$ and then taking trace on both sides, we obtain the equation for $`A(p^2)`$

$$A(p^2)=1+\frac{e^2}{4p^2}\frac{d^3q}{(2\pi )^3}tr[(i\gamma p)\gamma _\mu S(q)\mathrm{\Gamma }_\nu (p,q)D_{\mu \nu }(pq)].$$
(7)
In this paper, in order to satisfy the Ward-Takahashi identity, the vertex function  is taken to be

$$\mathrm{\Gamma }_\mu (p,q)=\gamma _\mu G(p^2,q^2).$$
(8)
The one-loop vacuum polarization  is

$$\mathrm{\Pi }(p^2)=\alpha /p$$
(9)
with $`\alpha =Ne^2/8.`$

Now, we can write the equation for $`B(p^2)`$ as:

$$B(p^2)=\frac{e^2}{4}\frac{d^3q}{(2\pi )^3}\frac{B(q^2)G(p^2,q^2)}{q^2A^2(q^2)+B^2(q^2)}tr[\gamma _\mu \gamma _\nu D_{\mu \nu }(pq)].$$
(10)
To look for DCSB in QED$`_3,`$ the nonlinear integral equations for $`A(p^2)`$ and $`B(p^2)`$ should be solved explicitly. It is easy to see that there is always a trivial solution $`B(p^2)0,`$ which means there is no mass generation of fermions and hence no DCSB. The conventional opinion claimed that a nontrivial solution $`B(p^2)`$ leads to DCSB. However, according to the well-known Goldstone theorem , the spectrum of physical particles must contain one massless spin-0 particle for each broken continuous symmetry. If the broken symmetry is a global chiral symmetry, the corresponding Goldstone boson is a fermion-antifermion ($`A`$ and $`\overline{A}`$) pseudoscalar bound state. When we study DCSB, the properties of the Goldstone boson must be taken into account. We have seen in Ref.  that the normalization condition of this bound state places a nontrivial constraint on the asymptotic form of $`A(p^2)`$ and $`B(p^2)`$ in quenched planar QED$`_4.`$ Now, we address this constraint in the three-dimensional QED.

This bound state is a nontrivial solution of the following BS equation in Euclidean space :

$$S_A^1(\frac{k}{2}+p)\chi _k(p)S_{\overline{A}}^1(\frac{k}{2}p)=e^2\frac{d^3q}{(2\pi )^3}\gamma _\mu \chi _k(q)\mathrm{\Gamma }_\nu (p,q)D_{\mu \nu }(pq).$$
(11)
Here, $`k/2=(p_A+p_{\overline{A}})/2`$ and $`p=(p_Ap_{\overline{A}})`$ are the center-of-mass 3-momentum and relative 3-momentum of the bound state, respectively. $`\chi _k(p)`$ is the bound state wave function expressed as a 4$`\times 4`$ matrix.

Combining Eq.(2) and Eq.(11), we have

$`\left[i\gamma \left({\displaystyle \frac{k}{2}}+p\right)A\left(\left({\displaystyle \frac{k}{2}}+p\right)^2\right)+B\left(\left({\displaystyle \frac{k}{2}}+p\right)^2\right)\right]`$ (12)
$`\times \chi _k(p)\left[i\gamma \left({\displaystyle \frac{k}{2}}p\right)A\left(\left({\displaystyle \frac{k}{2}}p\right)^2\right)B\left(\left({\displaystyle \frac{k}{2}}p\right)^2\right)\right]`$
$`=`$ $`e^2{\displaystyle \frac{d^3q}{(2\pi )^3}\gamma _\mu \chi _k(q)\mathrm{\Gamma }_\nu (p,q)D_{\mu \nu }(pq)}`$
Note that since we adopt four-component spinors, there are two matrices which anticommute with $`\gamma _0`$, $`\gamma _1`$ and $`\gamma _2`$. The mass term $`m\overline{\psi }\psi `$ breaks two chiral symmetries simultaneously. As a result there are two Goldstone bosons respectively coupling to $`\overline{\psi }\gamma _\mu \gamma _3\psi `$ and $`\overline{\psi }\gamma _\mu \gamma _5\psi `$. These Goldstone bound states both have the quantum number $`l^p=0^{}`$, which corresponds to the broken generator of the chiral group and vanishing 3-momentum. Their BS amplitudes are $`\chi _k(p)=\chi _{03}(p^2)\gamma _3`$ and $`\chi _k(p)=\chi _{05}(p^2)\gamma _5`$ respectively with $`k=0`$. It is easy to show that $`\chi _{03}(p^2)`$ and $`\chi _{05}(p^2)`$ satisfy the same BS equation and normalization condition, which give the same constraint on the nontrivial solution of the DS equation. For simplicity, we take $`\gamma _5`$ as the example and the result remains correct for $`\gamma _3`$.

Substituting $`\chi _k(p)=\chi _{05}(p^2)\gamma _5`$ and $`k=0`$ into Eq.(12), and taking trace on both sides, we finally obtain the BS equation for the Goldstone boson:

$$[p^2A^2(p^2)+B^2(p^2)]\chi _{05}(p^2)=\frac{e^2}{4}\frac{d^3q}{(2\pi )^3}\chi _{05}(q^2)G(p^2,q^2)tr[\gamma _\mu \gamma _\nu D_{\mu \nu }(pq)],$$
(13)
where $`A(p^2)`$ and $`B(p^2)`$ are solutions of Eq.(6) and Eq.(7) respectively. In the derivation of the right-hand side of this equation, we have used the approximation $`\mathrm{\Gamma }_\nu (p,q)=\gamma _\nu G(p^2,q^2).`$

The BS equation (13) is a homogeneous integral equation and hence can not determine the bound state completely, leaving an arbitrary multiplicative finite constant $`C.`$ To eliminate this uncertainty, an additional normalization condition of $`\chi _{05}(p^2)`$ is necessary. The search for such a normalization condition has a long history, and there are several approaches (different but equivalent) to normalize the bound state wave function . In this paper, for convenience, we shall utilize the form given by Suttorp . Rewriting the normalization condition in three dimension, we have

$`{\displaystyle }d^3ptr\{\overline{\chi }_k\left(p\right)[i\gamma ({\displaystyle \frac{k}{2}}+p)A\left(({\displaystyle \frac{k}{2}}+p)^2\right)+B\left(({\displaystyle \frac{k}{2}}+p)^2\right)]`$
$`\times \chi _k\left(p\right)[i\gamma ({\displaystyle \frac{k}{2}}p)A\left(({\displaystyle \frac{k}{2}}p)^2\right)B\left(({\displaystyle \frac{k}{2}}p)^2\right)]\}`$
$`=\lambda {\displaystyle \frac{dM_B}{d\lambda }}`$ (14)
where $`M_B`$ is the mass of the bound state and $`\lambda =e`$ is the coupling constant. Generally, it is sufficient to require that the integral at the left-hand side of Eq.(14) has a finite value . Setting $`k=0,`$ $`\chi _k(p)=\chi _{05}(p^2)\gamma _5`$ and using the identity $`\gamma _\mu \gamma _5=\gamma _5\gamma _\mu `$, the normalization condition can be written in the simple form

$$_0^{\mathrm{}}𝑑qq^2[q^2A^2(q^2)+B^2(q^2)]\chi _{05}^2(q^2)=finite.$$
(15)
Now we introduce the function

$$\varphi (p^2)=[p^2A^2(p^2)+B^2(p^2)]\chi _{05}(p^2),$$
(16)
then Eq.(13) becomes

$$\varphi (p^2)=\frac{e^2}{4}\frac{d^3q}{(2\pi )^3}\frac{\varphi (q^2)}{q^2A^2(q^2)+B^2(q^2)}G(p^2,q^2)tr[\gamma _\mu \gamma _\nu D_{\mu \nu }(pq)].$$
(17)
If there is a solution $`B(p^2)`$ of Eq.(10) (the equation for $`A(p^2)`$ always has a solution), then we also have a solution of Eq.(17) as $`\varphi (p^2)=cB(p^2),`$ with $`c`$ an arbitrary finite constant. According to Eq.(16), the Goldstone boson $`\chi _{05}(p^2)`$ is given by $`\varphi (p^2),`$ $`A(p^2)`$ and $`B(p^2)`$. Thus, we have constructed a relationship between two systems, one is described by DS equation and the other BS equation. Therefore, Eq.(16) and Eq.(17) may be regarded as another form of the Goldstone theorem.

Using Eq.(16), the normalization condition Eq.(15) can be written as

$$_0^{\mathrm{}}𝑑q\frac{q^2\varphi ^2(q^2)}{q^2A^2(q^2)+B^2(q^2)}=finite.$$
(18)
For chiral symmetry to be broken, the nontrivial solution of the DS equation (10) and the BS equation (13) must exist simultaneously. Besides, as a bound state wave function, each nontrivial solution of the BS equation satisfies a normalization condition. We can easily turn this condition into a constraint placed on the nontrivial solution of the DS equation, and now state a theorem.

Theorem: The necessary and sufficient condition for a nontrivial solution $`B(p^2)`$ to be chiral symmetry breaking solution is that it must satisfy, together with the solution $`A(p^2)`$, the condition

$$_0^{\mathrm{}}𝑑q\frac{q^2B^2(q^2)}{q^2A^2(q^2)+B^2(q^2)}=finite.$$
(19)
The proof of a similar theorem in the case of QED4 given in Ref.  is still valid in the current theory although the conditions in these two cases are different. Hence, we will not present the proof in this paper.

It is now straightforward to apply this result to the nontrivial solution of the DS equation. Although the approximations (Eq.(8) and Eq.(9)) that we take in this paper simplify the DS equation significantly, it appears impossible to obtain its completely analytical solution. To solve Eq.(10) analytically, further approximations are needed. Under the bare vertex and one loop vacuum-polarization approximations, Nash  considered the leading and next-to-leading terms in the 1/$`N`$ expansion of the kernels in the DS equation and showed that higher-order corrections do not alter the nature of the symmetry breaking. This result confirms the qualitative conclusion made earlier by Appelquist and coworkers  that there is a finite critical number of flavors $`N_c`$ above which the DS equation for $`B(p^2)`$ has no nontrivial solutions. Nash solved the DS equation for $`B(p^2)/A(p^2)`$ in low momentum region and showed that $`N_c`$ is gauge-independent. For $`N<N_c`$, considering only first order of 1/$`N`$ in this paper, $`B(p^2)/A(p^2)`$ has the form

$$\frac{B(p^2)}{A(p^2)}=p^{\frac{1}{2}}\mathrm{sin}\left(\frac{1}{2}\left[\frac{128}{3\pi ^2N}1\right]^{\frac{1}{2}}\left\{\mathrm{ln}\left[\frac{pA(0)}{B(0)}+\delta \right]\right\}\right),$$
(20)
in the infrared region, where $`\delta `$ is a phase and $`B(0)/A(0)`$ is finite. In the ultraviolet region ($`p\alpha `$), the mass function $`B(p^2)/A(p^2)`$ falls like 1/$`p^2`$ . This result was obtained via the operator-product expansion approach and does not depend on the 1/$`N`$ expansion .

In the current case, the condition (19) becomes

$$_0^{\mathrm{}}𝑑q\frac{q^2}{1+q^2A^2(q^2)/B^2(q^2)}=finite.$$
(21)
It is easy to see that there are no singularities in the integrand. The dangers those would cause Eq.(21) invalid come from the ultraviolet and infrared behaviors of $`B(p^2)/A(p^2)`$. Direct calculations show that the ultraviolet and infrared forms of $`B(p^2)/A(p^2)`$ satisfy this condition. Thus, at least to the lowest order in $`1/N`$ expansion, the nontrivial solution of the DS equation in QED3 does lead to DCSB and dynamical mass generation for fermions.

In conclusion, we have obtained a sufficient and necessary condition for a nontrivial solution of the DS equation to be chiral symmetry breaking solution. In the derivation of this condition, the Goldstone theorem and the normalization condition of the Goldstone bound state wave function play important roles. It is showed that the nontrivial solutions given by Nash and Appelquist $`et`$ $`al`$.  satisfy this condition. Therefore we see that QED in four dimension (in quenched planar approximation) and three dimension (under the approximations mentioned above) have different chiral phase structures. The former undergoes chiral phase transition only when a cutoff is present, but the latter exhibits DCSB in the continuum form when the number of fermions is less than a critical value. The origin of this difference is still unknown and subjected to the future studies.

It is claimed  that the intrinsic scale $`\alpha `$ divides QED3 into two parts and everything in high energy beyond $`\alpha `$ is rapidly damped. From this point of view, it seems reasonable to construct an equivalence between the continuous field theory, QED3, and the low-energy effective theory of antiferromagnet defined on discrete lattices. In some sense, $`\alpha `$ is the lattice constant of QED3 $`[8]`$.

One of us (G.L) would like to thank T.J. Du, H.J. Pan and T. Tu for helpful discussions. G.C. is supported by the National Science Foundation in China.