# Singularities and Pseudogaps in the Density of States of Peierls Chains

## Abstract

We develop a non-perturbative method to calculate the density of states (DOS) $`\rho (\omega )`$ of the fluctuating gap model describing the low-energy physics of electrons on a disordered Peierls chain. For a real order parameter field we calculate $`\rho (0)`$ (i.e. the DOS at the Fermi energy) exactly as a functional of the disorder for a chain of finite length $`L`$. Averaging $`\rho (0)`$ with respect to a Gaussian probability distribution of the fluctuating Peierls order parameter, we show that for $`L\mathrm{}`$ the average $`\rho (0)`$ diverges for any finite value of the correlation length above the Peierls transition. Pseudogap behavior emerges only if the Peierls order parameter is finite and sufficiently large.

PACS numbers: 71.23.-k, 02.50.Ey, 71.10.Pm

At low temperatures many quasi one-dimensional conductors become unstable and develop long-range charge-density-wave order, i.e. undergo a Peierls-transition. Within a mean field picture a finite value $`\mathrm{\Delta }_0`$ of the Peierls order parameter leads for frequencies $`|\omega |<|\mathrm{\Delta }_0|`$ to a gap in the electronic density of states (DOS) $`\rho (\omega )`$. To achieve a better understanding of the effect of fluctuations on the Peierls transition, Lee, Rice, and Anderson introduced the so-called fluctuating gap model (FGM). In this model the fluctuating part $`\stackrel{~}{\mathrm{\Delta }}(x)=\mathrm{\Delta }(x)\mathrm{\Delta }_0`$ of the order parameter is approximated by a Gaussian stochastic process with covariance

$$\stackrel{~}{\mathrm{\Delta }}(x)\stackrel{~}{\mathrm{\Delta }}(x^{})K(x,x^{})=K_0e^{|xx^{}|/\xi }.$$
(1)   
Here $`\mathrm{}`$ denotes averaging over the probability distribution of $`\mathrm{\Delta }(x)`$, $`K_0`$ is a positive constant, and $`\xi `$ is the order parameter correlation length. We assume that the field $`\mathrm{\Delta }(x)`$ is real, corresponding to a charge-density wave that is commensurate with the lattice.

Twenty years ago Sadovskii found an apparently exact algorithm to calculate the average DOS of the FGM. His calculations showed that for temperatures above the Peierls transition, in a regime where $`\xi `$ is large but finite, the average DOS exhibits a substantial suppression in the vicinity of the Fermi energy, a so-called pseudogap. The algorithm constructed by Sadovskii has also been applied in a different context to explain the weak pseudogap behavior in the underdoped cuprates. However, recently it has been pointed out that Sadovskii’s algorithm contains a subtle flaw and hence does not produce the exact DOS of the FGM. It is therefore important to compare this algorithm with limiting cases where $`\rho (\omega )`$ can be calculated without any approximation.

Besides the limit $`\xi \mathrm{}`$ where Sadovskii’s algorithm is indeed exact, there exists another non-trivial limit where the exact $`\rho (\omega )`$ is known: if in Eq.(1) we let $`\xi 0`$, $`K_0\mathrm{}`$, with $`K_0\xi D=\mathrm{const}`$, the right-hand side of Eq.(1) reduces to $`2D\delta (xx^{})`$. As shown by Ovchinnikov and Erikhman (OE), in the limit $`L\mathrm{}`$ (where $`L`$ is the length of the chain) the exact $`\rho (\omega )`$ can then be obtained from the stationary solution of a Fokker-Planck equation. For small $`\omega `$ and $`\mathrm{\Delta }_0=0`$ one finds $`\rho (\omega )|\omega \mathrm{ln}^3|\omega ||^1`$. Singularities of this type at the band center of a random Hamiltonian have been discovered by Dyson, and have recently also been found in one-dimensional spin-gap systems. It is important to note that in the FGM the singularity is a consequence of the charge conjugation symmetry of the underlying Dirac Hamiltonian, and is not related to concrete probability properties of $`\mathrm{\Delta }(x)`$. In particular, the singularity is not an artefact of the exactly solvable limit $`\xi 0`$ considered by OE. It is therefore reasonable to expect that for any $`\xi <\mathrm{}`$ the average DOS of the FGM exhibits a singularity at $`\omega =0`$. This general argument is in disagreement with Ref., where for large but finite $`\xi `$ a pseudogap (and hence no singularity) has been obtained. In this work we shall resolve this contradiction by calculating the average DOS at the Fermi energy (i.e. $`\rho (\omega =0)`$) exactly for arbitrary $`\xi `$.

The local DOS $`\rho (x,\omega )`$ of the FGM for a given realization of the disorder can be written as

$$\rho (x,\omega )=\pi ^1\mathrm{ImTr}[\sigma _3𝒢(x,x,\omega +i0)],$$
(2)   
where the $`2\times 2`$ matrix Green’s function $`𝒢`$ satisfies

$$\left[i_x+\omega \sigma _3i\mathrm{\Delta }(x)\sigma _2\right]𝒢(x,x^{},\omega )=\delta (xx^{})\sigma _0.$$
(3)   
Here $`\sigma _i`$ are the usual Pauli matrices and $`\sigma _0`$ is the $`2\times 2`$ unit matrix. Note that in Eq.(2) we have factored out a Pauli matrix $`\sigma _3`$, so that the differential operator $`i_x`$ in Eq.(3) is proportional to the unit matrix. To solve Eq.(3), we try the ansatz (suppressing for simplicity the frequency label)

$$𝒢(x,x^{})=U(x)𝒢_1(x,x^{})U^1(x^{}),$$
(4)   
where $`U(x)`$ is an invertible $`2\times 2`$ matrix. Eq.(4) resembles the transformation law of the comparator in non-Abelian gauge theories. In fact, Eq.(4) can be viewed as a gauge transformation which generalizes the Schwinger ansatz to the non-Abelian case. It is easy to show that the solution of Eq.(3) can indeed be written in the form (4) provided $`𝒢_1`$ and $`U`$ satisfy

$`\left[i_x+\omega \sigma _3\right]𝒢_1(x,x^{})`$ $`=`$ $`\delta (xx^{})\sigma _0,`$ (6)   
$`i_xU(x)=\omega [U(x)\sigma _3\sigma _3U(x)]+i\mathrm{\Delta }(x)\sigma _2U(x).`$   
Eq.(6) defines the Green’s function of free fermions, and can be solved trivially via Fourier transformation. The difficult part of the calculation is the solution of the matrix equation (6). We parameterize $`U(x)`$ as follows,

$$U(x)=e^{i\mathrm{\Phi }_+(x)\sigma _{}}e^{i\mathrm{\Phi }_{}(x)\sigma _+}e^{i\mathrm{\Phi }_3(x)\sigma _3},$$
(7)   
where $`\sigma _\pm =\frac{1}{2}[\sigma _1\pm i\sigma _2]`$, and the three functions $`\mathrm{\Phi }_\pm (x)`$, $`\mathrm{\Phi }_3(x)`$ have to be chosen such that $`U(x)`$ satisfies Eq.(6). A parameterization similar to Eq.(7) has recently been used by Schopohl to study the Eilenberger equations of superconductivity. We find that the ansatz (7) solves Eq.(6) if $`\mathrm{\Phi }_\pm (x)`$ and $`\mathrm{\Phi }_3(x)`$ satisfy

$`_x\mathrm{\Phi }_+`$ $`=`$ $`2i\omega \mathrm{\Phi }_++\mathrm{\Delta }(x)[1\mathrm{\Phi }_+^2],`$ (9)   
$`_x\mathrm{\Phi }_{}`$ $`=`$ $`2i\omega \mathrm{\Phi }_{}\mathrm{\Delta }(x)[12\mathrm{\Phi }_+\mathrm{\Phi }_{}],`$ (10)   
$`_x\mathrm{\Phi }_3`$ $`=`$ $`i\mathrm{\Delta }(x)\mathrm{\Phi }_+.`$ (11)   
Non-linear differential equations of the type (9) are called Riccati equations. The set of equations obtained by Schopohl has a similar structure but is not identical with Eqs.(9-11). Note that Eq.(9) involves only $`\mathrm{\Phi }_+`$. If we manage to obtain the solution $`\mathrm{\Phi }_+`$, Eqs.(10,11) become simple linear equations which can be solved exactly. From Eqs.(2,4) and (7) it is easy to see that the local DOS can be written as

$$\rho (x,\omega )=\pi ^1\mathrm{Re}R(x,\omega +i0),R=12\mathrm{\Phi }_+\mathrm{\Phi }_{}.$$
(12)   
Thus, to calculate the average DOS we have to average the product $`\mathrm{\Phi }_+\mathrm{\Phi }_{}`$ over the probability distribution of the field $`\mathrm{\Delta }(x)`$. In the limit where the right-hand-side of Eq.(1) reduces to $`2D\delta (xx^{})`$ we can use the fact that $`\mathrm{\Phi }_+`$ and $`\mathrm{\Phi }_{}`$ satisfy first order differential equations to express the average $`\mathrm{\Phi }_+\mathrm{\Phi }_{}`$ in terms of the solution of two coupled one-dimensional Fokker-Planck equations. Thus, our method leads to an algorithm for obtaining the exact $`\rho (\omega )`$ without using the node counting theorem. However, the Fokker-Planck equation obtained by OE within the phase formalism is easier to solve than our system of two coupled Fokker-Planck equations. Hence, for $`\delta `$-function correlated disorder our approach does not have any practical advantage.

On the other hand, if the disorder is not $`\delta `$-function correlated, probability distributions of physical quantities do in general not satisfy Fokker-Planck equations, and it is not so easy to perform controlled calculations or even obtain exact results. We now show that for real $`\mathrm{\Delta }(x)`$ the local DOS at the Fermi energy can be calculated exactly. To derive this result, let us introduce the complex vector

$$\stackrel{}{\psi }=\left(\begin{array}{c}\sqrt{2}(1\mathrm{\Phi }_+\mathrm{\Phi }_{})\mathrm{\Phi }_+\\ 12\mathrm{\Phi }_+\mathrm{\Phi }_{}\\ \sqrt{2}\mathrm{\Phi }_{}\end{array}\right)\left(\begin{array}{c}Z_+\\ R\\ Z_{}\end{array}\right).$$
(13)   
Note that by construction $`R^22Z_+Z_{}=1`$ for all $`x`$, and that according to Eq.(12) the second component of $`\stackrel{}{\psi }`$ is related to the local DOS. Using Eqs.(9,10) we find

$$_x\stackrel{}{\psi }=H(x)\stackrel{}{\psi },H(x)=2i\omega J_3+2\mathrm{\Delta }(x)J_1,$$
(14)   
where $`J_i`$ are spin $`J=1`$ operators in the representation

$$J_3=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right),J_1=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right).$$
(15)   
Eq.(14) is a linear multiplicative stochastic differential equation. Formally this equation looks like the imaginary time Schrödinger equation for a $`J=1`$ quantum spin in a random magnetic field, with $`x`$ playing the role of imaginary time. Although the operator $`H`$ in Eq.(14) is not hermitian, we may perform an analytic continuation to imaginary frequencies ($`\omega =iE`$) to obtain a hermitian spin Hamiltonian. We thus arrive at the remarkable result that the average DOS of the FGM can be obtained from the average state-vector of a $`J=1`$ spin in a random magnetic field. Our non-linear transformation (13) is well known in the quantum theory of magnetism: with the formal identification $`RJ_3`$, $`\sqrt{2}Z_\pm J_\pm `$, $`\sqrt{2}\mathrm{\Phi }_{}b^{}`$, and $`\sqrt{2}\mathrm{\Phi }_+b`$, Eq.(13) is precisely the Dyson-Maleev transformation, which expresses the spin operators in terms of boson operators $`b`$,$`b^{}`$.

The solution of Eq.(14) with initial condition $`\stackrel{}{\psi }(0)=\stackrel{}{\psi }_0`$ is $`\stackrel{}{\psi }(x)=S(x)\stackrel{}{\psi }_0`$, where the $`S`$-matrix is

$$S(x)=𝒯\mathrm{exp}\left[_0^x𝑑x^{}H(x^{})\right].$$
(16)   
Here $`𝒯`$ is the usual time ordering operator. The proper choice of boundary conditions requires some care. For simplicity, let us assume that $`\mathrm{\Delta }(x)`$ is non-zero only in a finite interval $`0xL`$. Outside this regime we find that $`Z_\pm (x)=\mathrm{exp}[2i(\omega +i0)(xx_0)]Z_\pm (x_0)`$ is the solution of Eq.(14). Physically it is clear that exponentially growing solutions are forbidden, which requires $`Z_{}(0)=Z_+(L)=0`$ and implies $`R=1`$ for $`x0`$ and $`xL`$. We conclude that at the boundaries

$$\stackrel{}{\psi }_0=\left(\begin{array}{c}Z_+(0)\\ 1\\ 0\end{array}\right),\stackrel{}{\psi }(L)=\left(\begin{array}{c}0\\ 1\\ Z_{}(L)\end{array}\right),$$
(17)   
where $`Z_+(0)`$ and $`Z_{}(L)`$ are determined by $`\stackrel{}{\psi }(L)=S(L)\stackrel{}{\psi }_0`$. This implies $`Z_+(0)=S_{12}(L)/S_{11}(L)`$. Because the matrix elements of $`S(L)`$ depend on the disorder, the initial vector $`\stackrel{}{\psi }_0`$ is stochastic. Note that in textbook discussions of multiplicative stochastic differential equations one often assumes deterministic initial conditions. After simple algebra we obtain for the second component of $`\stackrel{}{\psi }(x)=S(x)\stackrel{}{\psi }_0`$ in the interval $`0xL`$

$$R(x)=S_{22}(x)S_{21}(x)S_{12}(L)/S_{11}(L).$$
(18)   
In general we have to rely on approximations to calculate the time-ordered exponential in Eq.(16). However, there are two special cases where $`S(x)`$ can be calculated exactly. The first is obvious: if $`\mathrm{\Delta }(x)=\mathrm{\Delta }_0`$ is independent of $`x`$, our spin Hamiltonian $`H(x)`$ is constant, so that the time-ordering operator is not necessary. Eq.(18) can then be evaluated exactly for arbitrary $`L`$. If we take the limit $`L\mathrm{}`$ holding $`x/L`$ fixed, we recover the well known square root singularity at the band edges,

$$\underset{L\mathrm{}}{lim}\rho (x,\omega )=\frac{\mathrm{\Theta }(\omega ^2\mathrm{\Delta }_0^2)|\omega |}{\pi \sqrt{\omega ^2\mathrm{\Delta }_0^2}},\mathrm{\hspace{0.33em}\hspace{0.33em}0}<x/L<1.$$
(19)   
There exists another, more interesting limit where $`S(x)`$ can be calculated exactly. Obviously, at $`\omega =0`$ the direction of the magnetic field in our spin Hamiltonian (14) is constant. Although in this case $`H(x)`$ is $`x`$-dependent, we may omit the time-ordering operator in Eq.(16). After straightforward algebra we obtain from Eq.(18)

$$R(x)=\mathrm{cosh}[A(x)B(x)]/\mathrm{cosh}[A(x)+B(x)],$$
(20)   
where $`A(x)=_0^x𝑑x^{}\mathrm{\Delta }(x^{})`$ and $`B(x)=_x^L𝑑x^{}\mathrm{\Delta }(x^{})`$. In Eq.(20) it is understood that $`R(x)`$ stands for $`R(x,i0)`$, so that $`\rho (x,0)=\pi ^1R(x)`$. We have thus succeeded to calculate the local DOS $`\rho (x,\omega =0)`$ of the FGM at the Fermi energy for a given realization of the disorder. The special symmetries of random Dirac fermions at $`\omega =0`$ have recently been used by Shelton and Tsvelik to calculate the statistics of the corresponding wave-functions.

To calculate the disorder average of Eq.(20), we introduce $`P(x;a,b)=\delta (aA(x))\delta (bB(x))`$ and write 

$$R(x)=_{\mathrm{}}^{\mathrm{}}𝑑a_{\mathrm{}}^{\mathrm{}}𝑑bP(x;a,b)\frac{\mathrm{cosh}(ab)}{\mathrm{cosh}(a+b)}.$$
(21)   
Assuming that the probability distribution of $`\mathrm{\Delta }(x)`$ is Gaussian with average $`\mathrm{\Delta }_0`$ and covariance $`K(x,x^{})`$, the joint distribution $`P(x;a,b)`$ can be calculated exactly. The integration over the difference $`ab`$ in Eq.(21) can then be performed, and we obtain for the average local DOS $`\rho (x,0)=\pi ^1R(x)`$,

$`\rho (x,0)`$ $`=`$ $`{\displaystyle \frac{e^{\alpha (x)}}{\pi \sqrt{2\pi \sigma ^2}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑s\mathrm{exp}\left[s^2/(2\sigma ^2)\right]`$ (22)   
$`\times `$ $`{\displaystyle \frac{\mathrm{cosh}[\beta (x)s+\mathrm{\Delta }_0(2xL)]}{\mathrm{cosh}[s+\mathrm{\Delta }_0L]}}.`$ (23)   
Here $`\sigma ^2=C_1(L)`$, $`\alpha (x)=2[C_1(x)C_2(x)C_3^2(x)]/\sigma ^2`$, and $`\beta (x)=[C_1(x)C_2(x)]/\sigma ^2`$, with

$$C_1(x)=_0^x𝑑x^{}_0^x𝑑x^{\prime \prime }K(x^{},x^{\prime \prime }).$$
(24)   
$`C_2(x)`$ is defined by replacing the range of the integrals in Eq.(24) by the interval $`[x,L]`$, and $`C_3(x)`$ is obtained by choosing the interval $`[0,x]`$ for the $`x^{}`$\- and $`[x,L]`$ for the $`x^{\prime \prime }`$-integration.

We now specify $`K(x,x^{})`$ to be of the form (1). Then $`\sigma ^2`$, $`\alpha (x)`$, and $`\beta (x)`$ are easily calculated. It is convenient to introduce the dimensionless parameters $`\stackrel{~}{L}=2K_0\xi L`$, $`\lambda =\xi /L`$, and $`\nu =\mathrm{\Delta }_0/(2K_0\xi )`$. Numerical results for $`\rho (x,0)`$ are shown in Fig.1. Due to symmetry with respect to $`x=L/2`$ the local DOS assumes an extremum at $`x=L/2`$, which in the limit $`L\mathrm{}`$ approaches either zero or infinity. Using the fact that at $`x=L/2`$ the $`\mathrm{cosh}`$ in the numerator of Eq.(23) is unity, we obtain

$$\rho (L/2,0)=\frac{e^{\frac{\stackrel{~}{L}}{2}f(\lambda )}}{\pi \sqrt{2\pi \sigma ^2}}_{\mathrm{}}^{\mathrm{}}𝑑s\frac{\mathrm{exp}\left[s^2/(2\sigma ^2)\right]}{\mathrm{cosh}[s+\nu \stackrel{~}{L}]},$$
(25)   
where $`\sigma ^2=\stackrel{~}{L}[1\lambda (1e^{1/\lambda })]`$ and

$$f(\lambda )=1\lambda [34e^{1/(2\lambda )}+e^{1/\lambda }].$$
(26)   
As shown in Fig.2, $`f(\lambda )`$ is positive and monotonically decreasing. Let us first consider the case $`\nu =0`$. This corresponds to the model discussed by Sadovskii with real $`\mathrm{\Delta }(x)`$. For $`\stackrel{~}{L}1`$ the $`s`$-integration in Eq.(25) is easily done, and we find $`\rho (\frac{L}{2},0)\stackrel{~}{L}^{1/2}\mathrm{exp}[\frac{\stackrel{~}{L}}{2}f(\lambda )]`$. Keeping in mind that for any finite $`\xi `$ the parameter $`\lambda =\xi /L`$ vanishes for $`L\mathrm{}`$, it is obvious that in this limit $`\rho (L/2,0)`$ is infinite. From Eq.(23) it is easy to show numerically that this is also true for $`lim_L\mathrm{}\rho (x,0)`$ in the open interval $`0<x/L<1`$. We have thus proven that for $`L\mathrm{}`$ and arbitrary $`\xi <\mathrm{}`$ the average DOS of the FGM is infinite at the Fermi energy, in agreement with general symmetry arguments.

For finite $`\nu `$ a careful analysis of Eq.(25) shows that there exists a critical value $`\nu _c(\lambda )`$ such that for $`|\nu |>\nu _c`$ the local DOS $`\rho (L/2,0)`$ scales to zero in the thermodynamic limit, and a pseudogap emerges. We obtain

$$\nu _c(\lambda )=[1\lambda (1e^{1/\lambda })]^{1/2}[f(\lambda )]^{1/2},$$
(27)   
see Fig.2. For $`\lambda =0`$ Eq.(27) yields $`\nu _c(0)=1`$. At the first sight this seems to contradict the result of OE, who found pseudogap behavior already for $`|\nu |>1/2`$. One should keep in mind, however, that we have set $`\omega =0`$ before taking the limit $`L\mathrm{}`$, while in Ref. these limits are taken in the opposite order. The non-commutativity of these limits is well known from the calculation of the local DOS of the Tomonaga-Luttinger model with a boundary. Interestingly, for $`K(x,x^{})=2D\delta (xx^{})`$ there exists a regime $`1/2<\nu <1`$ where in the thermodynamic limit $`\rho (\omega )`$ is discontinuous at $`\omega =0`$. This follows from the fact that according to OE $`\rho (\omega )|\omega |^{2\nu 1}0`$ for $`\omega 0`$, whereas we have shown that $`\rho (0)=\mathrm{}`$. For $`\nu =0`$ and large but finite $`\xi `$ we conjecture that the frequency-dependence of $`lim_L\mathrm{}\rho (\omega )`$ is qualitatively similar to the behavior found by Fabrizio and Mélin for random Dirac fermions with a special type of disorder. Specifically, we expect that for frequencies exceeding a certain crossover frequency $`\omega ^{}`$ the average DOS of the FGM shows pseudogap behavior, which is correctly predicted by Sadovskii’s algorithm. However, this algorithm misses the Dyson singularity, which emerges for frequencies $`|\omega |\stackrel{<}{}\omega ^{}`$ for any finite value of $`\xi `$.

In summary, we have developed a non-perturbative method to calculate the Green’s function of the FGM. Our main result is the proof of the existence of a singularity in $`\rho (0)`$ for any finite value of the correlation length $`\xi `$ as long as $`\mathrm{\Delta }(x)`$ is real and $`|\mathrm{\Delta }_0|`$ is sufficiently small. For $`\mathrm{\Delta }_0=0`$ we have shown that with open boundary conditions $`\rho (\frac{L}{2},0)\mathrm{exp}[K_0\xi Lf(\xi /L)]`$, where $`f(0)=1`$. Moreover, if we let $`L\mathrm{}`$ keeping $`x/L(0,1)`$ fixed, $`\rho (x,0)`$ exhibits a similar singularity. Thus, disorder pushes states from high energies to the band center. For finite $`\mathrm{\Delta }_0`$ this effect competes with the suppression of the DOS due to long-range order. In the incommensurate case (where $`\mathrm{\Delta }(x)`$ is complex and in Eq. (3) we should replace $`i\sigma _2\mathrm{\Delta }\sigma _+\mathrm{\Delta }\sigma _{}\mathrm{\Delta }^{}`$ ) it is known that $`\rho (0)`$ is finite in the white-noise limit. We expect that this remains true for arbitrary $`\xi `$. In fact, there exists numerical evidence that in the incommensurate case $`\rho (\omega )`$ can be accurately calculated from Sadovskii’s algorithm. For a comparison with experiments one should keep in mind that any violation of the perfect charge conjugation symmetry will wash out the singularity at $`\omega =0`$. It is therefore unlikely that the singularity is visible in realistic materials, although an enhancement might survive. The fact that the singular behavior of $`\rho (0)`$ in the FGM with real $`\mathrm{\Delta }(x)`$ is only destroyed if $`|\nu |=|\mathrm{\Delta }_0|/(2K_0\xi )`$ exceeds a finite critical value implies that in commensurate Peierls chains true pseudogap behavior should emerge gradually below the Peierls transition when the order parameter $`\mathrm{\Delta }_0`$ is sufficiently large. Lee, Rice, and Anderson came to a similar conclusion within perturbation theory.

We thank M. V. Sadovskii, K. Schönhammer, R. Hayn, R. H. McKenzie, and H. Monien for discussions and comments. This work was financially supported by the DFG (Grants No. Ko 1442/3-1 and Ko 1442/4-1).