text:
# Dynamical density-matrix renormalization-group method

## I Introduction

The density-matrix renormalization group (DMRG) steve ; dmrgbook  is a very successful numerical methods for calculating static properties of ground states and low-lying eigenstates in quantum many-body systems. For low-dimensional strongly correlated systems DMRG is as accurate as exact diagonalization techniques but can be used to study much larger systems than with those techniques (currently, up to $`10^3`$ sites). Using a finite-size scaling analysis it is thus possible to determine the static properties of a system in the thermodynamic limit with great accuracy.

The calculation of dynamical properties and higher energy excitations with DMRG has proved to be more difficult. Several approaches have been proposed but calculations have been carried out successfully for few problems only. The simplest of these methods is the Lanczos DMRG. karen ; till  In practice, this method gives accurate results for the first few moments of a dynamical spectrum. Therefore, it works well for simple discrete spectra made of a few peaks but it usually fails for more complicated spectra. An alternate method for calculating dynamical properties is the correction vector DMRG. pati ; till  Contrary to the Lanczos DMRG, this method can describe complex or dense spectra accurately. Nevertheless, there have been few applications till2 ; brune  of correction vector DMRG until now because this method is more difficult and requires significantly more computer resources than the usual DMRG method for calculating static properties at low energy.

In this paper I describe a simple and efficient method, called the dynamical DMRG (DDMRG), for calculating dynamical properties and excited states with DMRG. This approach is based on a variational principle for dynamical correlation functions and the related excited states. The variational principle is essentially an elegant formulation of the correction vector technique. Because of the variational formulation, however, the DDMRG method is easier to use and significantly more accurate than the correction vector DMRG method.

While the spectrum of a finite system is necessarily discrete, continuous excitation bands are often found in the thermodynamic limit. It is possible to broaden finite-size spectra to simulate the continuum of an infinite-system spectrum. Usually, the broadening is arbitrarily large and no systematic quantitative analysis of finite-size effects is performed for the spectrum. Here I show that dynamical properties of infinite systems can be obtained reliably using an appropriate finite-size scaling analysis. The key to the analysis is the use of a broadening which scales systematically with the system size.

The DDMRG method and the finite-size scaling technique for dynamical spectrum have already been successfully used to investigate the optical properties of one-dimensional Mott insulators. eric ; fabian  Here I apply these techniques to the calculation of the optical conductivity in the one-dimensional Peierls-Hubbard model of conjugated polymers. dionys  Much effort has been devoted to understanding the optical properties of these materials. In particular, the optical conductivity of the Peierls-Hubbard model has been studied extensively and analytical results have been obtained for various special cases, such as the Mott-Hubbard insulator and the Peierls (band) insulator limits. Leaving out these special limits the Peierls-Hubbard model describes a Mott-Peierls insulator and its optical properties are still poorly understood. In this paper I show that DDMRG can accurately reproduce the known analytical results for the optical spectrum in the thermodynamic limit. Then I investigate the optical conductivity of a Mott-Peierls insulator using DDMRG and show that it displays specific features.

The paper is organized as follows: The variational principle for dynamical correlation functions and related excited states is presented in the next section. I describe the dynamical DMRG method in Sec. III. The finite-size scaling analysis is presented in Sec. IV. I report and discuss the results for the optical conductivity of the Peierls-Hubbard model in Sec. V. Finally, I conclude in the last section.

## II Variational principle

The dynamic response of a quantum system to a time-dependent perturbation is often given by dynamical correlation functions such as

$$G_A(\omega +i\eta )=\frac{1}{\pi }\psi _0|A^{}\frac{1}{E_0+\omega +i\eta H}A|\psi _0,$$
(1)
where $`H`$ is the time-independent Hamiltonian of the system, $`E_0`$ and $`|\psi _0`$ are its ground-state energy and wavefunction, $`A`$ is the quantum operator corresponding to the physical quantity which is analyzed, and $`A^{}`$ is the Hermitian conjugate of $`A`$. A small real number $`\eta >0`$ is used in the calculation to shift the poles of the correlation function into the complex plane. (I set $`\mathrm{}=1`$ in sections II to IV.)

In general, we are interested in calculating the imaginary part of the correlation function

$`I_A(\omega +i\eta )`$ $`=`$ $`\text{Im}G_A(\omega +i\eta )`$
$`=`$ $`{\displaystyle \frac{1}{\pi }}\psi _0|A^{}{\displaystyle \frac{\eta }{(E_0+\omega H)^2+\eta ^2}}A|\psi _0`$
in the limit $`\eta 0`$

$$I_A(\omega )=\underset{\eta 0}{lim}I_A(\omega +i\eta )0.$$
(3)
It should be noted that the spectrum $`I_A(\omega +i\eta )`$ for any finite $`\eta >0`$ can be calculated from the spectrum $`I_A(\omega )`$ by convolution with a Lorentzian distribution

$$I_A(\omega +i\eta )=C_\eta [I_A(\omega )]>0,$$
(4)
where I use the notation $`C_\eta [f(\omega )]`$ to represent the convolution of a spectral function $`f(\omega )`$ with a Lorentzian distribution of width $`\eta `$

$$C_\eta [f(\omega )]=_{\mathrm{}}^+\mathrm{}𝑑\omega ^{}f(\omega ^{})\frac{1}{\pi }\frac{\eta }{(\omega \omega ^{})^2+\eta ^2}.$$
(5)
The moments of the spectrum $`I_A(\omega )`$ fulfill sum rules such as

$`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega I_A(\omega )`$ $`=`$ $`\psi _0|A^{}A|\psi _0,`$
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega I_A(\omega )\omega `$ $`=`$ $`\psi _0|A^{}[H,A]|\psi _0,`$ (6)
$`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑\omega I_A(\omega )\omega ^2`$ $`=`$ $`\psi _0|[A^{},H][H,A]|\psi _0,`$
where $`[A,B]=ABBA`$.

A dynamical correlation function (1) can be calculated using the correction vector method. The correction vector associated with $`G_A(\omega +i\eta )`$ is defined by

$$|\psi _A(\omega +i\eta )=\frac{1}{E_0+\omega +i\eta H}|A,$$
(7)
where $`|A=A|\psi _0`$. If the correction vector is known, the dynamical correlation function can be calculated directly

$$G_A(\omega +i\eta )=\frac{1}{\pi }A|\psi _A(\omega +i\eta ).$$
(8)
To calculate a correction vector one first solves an inhomogeneous linear equation

$$\left[(E_0+\omega H)^2+\eta ^2\right]|\psi =\eta |A,$$
(9)
which always has a unique solution $`|\psi =|Y_A(\omega +i\eta )`$ for $`\eta 0`$. The correction vector is then given by

$$|\psi _A(\omega +i\eta )=|X_A(\omega +i\eta )+i|Y_A(\omega +i\eta )$$
(10)
with

$$|X_A(\omega +i\eta )=\frac{HE_0\omega }{\eta }|Y_A(\omega +i\eta ).$$
(11)
One should note that the states $`|X_A(\omega +i\eta )`$ and $`|Y_A(\omega +i\eta )`$ are complex if the state $`|A`$ is not real, but they always determine the real part and imaginary part of the dynamical correlation function $`G_A(\omega +i\eta )`$, respectively,

$`\text{Re}G_A(\omega +i\eta )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}A|X_A(\omega +i\eta ),`$ (12a)
$`\text{Im}G_A(\omega +i\eta )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}A|Y_A(\omega +i\eta ).`$ (12b)
The derivatives of the real and imaginary parts can also be calculated from these states

$`{\displaystyle \frac{d}{d\omega }}\text{Re}G_A(\omega +i\eta )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}[X_A(\omega +i\eta )|X_A(\omega +i\eta )`$ (13)
$`Y_A(\omega +i\eta )|Y_A(\omega +i\eta )]`$
$`{\displaystyle \frac{d}{d\omega }}\text{Im}G_A(\omega +i\eta )`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}X_A(\omega +i\eta )|Y_A(\omega +i\eta ).`$
A well-established approach for solving an inhomogeneous linear equation (9) is to formulate it as a minimization problem. One considers the functional

$`W_{A,\eta }(\omega ,\psi )`$ $`=`$ $`\psi |(E_0+\omega H)^2+\eta ^2|\psi `$ (14)
$`+\eta A|\psi +\eta \psi |A.`$
For any $`\eta 0`$ and a fixed frequency $`\omega `$ this functional has a well-defined and non-degenerate minimum for the quantum state which is solution of Eq. (9)

$$|\psi _{\text{min}}=|Y_A(\omega +i\eta ).$$
(15)
It is easy to show that the value of the minimum is related to the imaginary part of the dynamical correlation function

$$W_{A,\eta }(\omega ,\psi _{\text{min}})=\pi \eta I_A(\omega +i\eta ).$$
(16)
Therefore, the calculation of spectral functions can be formulated as a minimization problem. To determine $`I_A(\omega +i\eta )`$ at any frequency $`\omega `$ and for any $`\eta >0`$, one minimizes the corresponding functional $`W_{A,\eta }(\omega ,\psi )`$. Once this minimization has been carried out, the real part of the correlation function $`G_A(\omega +i\eta )`$ can be calculated using Eqs. (11) and (12a) if necessary. This is the variational principle for dynamical correlation functions.

It is clear that if we can calculate $`|Y_A(\omega +i\eta )`$ exactly, this variational formulation is completely equivalent to the correction vector method. However, if we can only calculate an approximate solution with an error of the order $`ϵ1`$, $`|\psi =|Y_A(\omega +i\eta )+ϵ|\varphi `$ with $`\varphi |\varphi =1`$, the variational formulation is more accurate. In the correction vector method the error in the spectrum $`I_A(\omega +i\eta )`$ calculated with Eq. (12b) is also of the order of $`ϵ`$. In the variational approach it is easy to show that the error in the value of the minimum $`W_{A,\eta }(\omega ,\psi _{\text{min}})`$, and thus in $`I_A(\omega +i\eta )`$, is of the order of $`ϵ^2`$. With both methods the error in the real part of $`G_A(\omega +i\eta )`$ is of the order of $`ϵ`$.

One can write the function $`I_A(\omega )`$ in the spectral form (or Lehmann representation)

$$I_A(\omega )=\underset{n}{}|\psi _n|A|\psi _0|^2\delta (\omega +E_0E_n),$$
(17)
where $`|\psi _0`$ is the ground state, $`|\psi _n,n>1,`$ denotes the other eigenstates of $`H`$, and $`E_0,E_n`$ are their respective eigenenergies. Obviously, only the eigenstates with a finite matrix element $`\psi _n|A|\psi _00`$ contribute to the spectrum and here we are only interested in those excited states. In the correction vector method the excitation energies $`E_nE_0`$ and the spectral weights $`|\psi _n|A|\psi _0|^2`$ can be obtained from the poles of $`G_A(\omega +i\eta )`$. The corresponding wavefunctions $`|\psi _n`$ can be calculated by taking the $`\eta 0`$ limit of the correction vectors

$$|\psi _n\underset{\eta 0}{lim}|Y_A(E_nE_0+i\eta ).$$
(18)
The excited states contributing to $`G_A(\omega +i\eta )`$ correspond to the local maxima of the spectrum $`I_A(\omega +i\eta )`$ for small enough $`\eta >0`$. Therefore, they can also be obtained by minimization of the functional $`W_{A,\eta }(\omega ,\psi )`$ with respect to both $`\omega `$ and $`\psi `$. The local minima of $`W_{A,\eta }(\omega ,\psi )`$ are given by the conditions

$`\omega _{\text{min}}+E_0`$ $`=`$ $`{\displaystyle \frac{\psi _{\text{min}}|H|\psi _{\text{min}}}{\psi _{\text{min}}|\psi _{\text{min}}}},`$
$`|\psi _{\text{min}}`$ $`=`$ $`|Y_A(\omega _{\text{min}}+i\eta ).`$ (19)
In the limit $`\eta 0`$, $`\omega _{\text{min}}+E_0`$ tends to the energy $`E_n`$ of an eigenstate with finite spectral weight, $`|\psi _{\text{min}}`$ is equal to the corresponding eigenstate $`|\psi _n`$ up to a normalization constant, and $`W_{A,\eta }(\omega _{\text{min}},\psi _{\text{min}})`$ tends to the spectral weight $`|\psi _n|A|\psi _0|^2`$. This is the variational principle for excited states contributing to a dynamical correlation function $`G_A(\omega +i\eta )`$.

Again this variational formulation is completely equivalent to the correction vector method if $`|Y_A(\omega +i\eta )`$ can be calculated exactly. In an approximate calculation, however, errors in the eigenenergies and spectral weights are of the order of $`ϵ`$ with the correction vector method, while they are of the order of $`ϵ^2`$ with the variational formulation.

## III Dynamical DMRG

DMRG is a numerical method for calculating the properties of lattice quantum many-body systems. It is described in detail in several publications (for instance, see Refs. steve,  and dmrgbook, ). DMRG can be considered as a variational approach. The system energy

$$E(\psi )=\frac{\psi |H|\psi }{\psi |\psi }$$
(20)
is minimized in a variational subspace (the DMRG basis) of the system Hilbert space to find the ground-state wavefunction $`|\psi _0`$ and energy $`E_0=E(\psi _0)`$. If the ground-state wavefunction is calculated with an error of the order of $`ϵ1`$ (i.e., $`|\psi =|\psi _0+ϵ|\varphi `$ with $`\varphi |\varphi =1`$), the energy obtained is an upper bound to the exact result and the error in the energy is of the order of $`ϵ^2`$ (as in all variational approaches). In principle, the DMRG energy error is proportional to the weight of the density-matrix eigenstates discarded in the renormalization procedure. This discarded weight can be reduced by increasing the number $`m`$ of density-matrix eigenstates kept in the calculation, which corresponds to an increase of the variational subspace dimension. Therefore, the energy error systematically decreases with increasing $`m`$ in a DMRG calculation.

The DMRG procedure used to minimize the energy functional (20) can also be used to minimize the functional $`W_{A,\eta }(\omega ,\psi )`$ and thus to calculate the dynamical correlation function $`G_A(\omega +i\eta )`$. I call this approach the dynamical DMRG method. The minimization of the functional $`W_{A,\eta }(\omega ,\psi )`$ is easily integrated into the usual DMRG algorithm. At every step of a DMRG sweep through the system lattice, a superblock representing the system is built and the following calculations are performed in the the superblock subspace:

1. The energy functional $`E(\psi )`$ is minimized using a standard iterative algorithm for the eigenvalue problem. This yields the ground state vector $`|\psi _0`$ and its energy $`E_0`$ in the superblock subspace.
2. The state $`|A`$ is calculated.
3. The functional $`W_{A,\eta }(\omega ,\psi )`$ is minimized using an iterative minimization algorithm. This gives the first part of the correction vector $`|Y_A(\omega +i\eta )`$ and the imaginary part $`I_A(\omega +i\eta )`$ of the dynamical correlation function.
4. The second part $`|X_A(\omega +i\eta )`$ of the correction vector is calculated using Eq. (11).
5. The real part and the derivatives of the dynamical correlation function can be calculated from Eqs. (12a) and (13), respectively.
6. The four states $`|\psi _0`$, $`|A`$, $`|Y_A(\omega +i\eta )`$, and $`|X_A(\omega +i\eta )`$ are included as target in the density-matrix renormalization to build a new superblock at the next step.

The robust finite-system DMRG algorithm must be used to perform several sweeps through a lattice of fixed size. Sweeps are repeated until the procedure has converged to the minimum of both functionals $`E(\psi )`$ and $`W_{A,\eta }(\omega ,\psi )`$.

To obtain the dynamical correlation function $`G_A(\omega +i\eta )`$ over a range of frequencies, one has to repeat this calculation for several frequencies $`\omega `$. If the DDMRG calculations are performed independently, the computational effort is roughly proportional to the number of frequencies. It is also possible to carry out a DDMRG calculation for several frequencies simultaneously, including several states $`|X_A(\omega +i\eta )`$ and $`|Y_A(\omega +i\eta )`$ with different frequencies $`\omega `$ as target. The optimal number of different frequencies to be included in a single calculation depends strongly on the problem studied and the computer used. As calculations for different frequencies are essentially independent, it would be easy and very efficient to perform these calculations on a parallel computer.

If one performs a DDMRG calculation for two close frequencies $`\omega _1`$ and $`\omega _2`$ simultaneously, it is possible to calculate the dynamical correlation function for additional frequencies $`\omega `$ between $`\omega _1`$ and $`\omega _2`$ without including the corresponding states $`|X_A(\omega +i\eta )`$ and $`|Y_A(\omega +i\eta )`$ as target in the density-matrix renormalization. This approach can significantly reduce the computer time necessary to determine the spectrum over a frequency range but the results obtained for $`\omega \omega _1,\omega _2`$ are less accurate and not always reliable, as the DMRG basis is optimized for the frequencies $`\omega _1`$ and $`\omega _2`$ only. (A similar technique is the calculation of spectra with the Lanczos algorithm in the DMRG basis optimized for a pair of correction vectors, see Ref. till, .) Alternatively, between the frequencies for which $`G_A(\omega +i\eta )`$ is determined directly with DDMRG, we can calculate the dynamical correlation function by interpolation using the DDMRG data for the function and its derivative.

If a complete spectrum $`I_A(\omega +i\eta )`$ has been obtained, it is possible to calculate the moments of the spectral distribution \[the left-hand-side of Eq. (6)\]. The first few moments \[the right-hand-side of Eq. (6)\] can be calculated accurately using the Lanczos DMRG method. karen ; till  This provides an independent check of DDMRG results. \[Note that only the first sum rule (6) is satisfied exactly for $`\eta >0`$.\]

To calculate individual excited states contributing to the spectrum in a given frequency range $`(\omega _1,\omega _2`$), one includes a minimization of $`W_{A,\eta }(\omega ,\psi )`$ with respect to $`\omega `$ ($`\omega _1<\omega <\omega _2`$) in the third step of the DDMRG algorithm described above. In this case, $`|Y_A(\omega _{\text{min}}+i\eta )`$ and $`|X_A(\omega _{\text{min}}+i\eta )`$ are included as target in the sixth step. The parameter $`\eta `$ must be much smaller than the distance $`E_{n+1}E_n`$ between two successive eigenstates contributing to the dynamical correlation function or must decrease during the calculation until the desired accuracy is obtained. To make the procedure robust it is necessary to simultaneously target a second correction vector $`|\psi _A(\omega +i\eta )`$ with a fixed frequency and a parameter $`\eta `$ of the order of the frequency range. Typically, I use $`\omega =(\omega _1+\omega _2)/2`$ and $`\eta =(\omega _2\omega _1)/4`$.

Because of the variational principle one naively expects that the DDMRG results for $`I_A(\omega +i\eta )`$ must converge monotonically from below to the exact result as the number $`m`$ of density-matrix eigenstates is increased. In practice, the convergence is not always regular because of two approximations made to calculate the functional $`W_{A,\eta }(\omega ,\psi )`$ in a DMRG basis. First, the ground-state energy $`E_0`$ and the state $`|A`$ used in the definition (14) of $`W_{A,\eta }(\omega ,\psi )`$ are not known exactly but calculated with DMRG. If the number $`m`$ of density matrix eigenstates is increased, $`E_0`$ and $`|A`$ are modified (they become progressively more accurate) and the functional $`W_{A,\eta }(\omega ,\psi )`$ is changed, which can affect its minimum arbitrarily. We also note that errors of the order of $`ϵ`$ in $`E_0`$ or $`|A`$ result in errors of the same order in $`I_A(\omega +i\eta )`$. Therefore, to observe a regular convergence with increasing $`m`$ and to obtain accurate results for $`I_A(\omega +i\eta )`$, it is necessary in the first place to determine the ground state and the state $`|A`$ with great precision (and thus to include the state $`|A`$ as target).

To calculate the functional $`W_{A,\eta }(\omega ,\psi )`$ in the third step of the DDMRG algorithm, one needs an effective representation of the operator $`(HE_0\omega )^2`$ in the superblock subspace

$$[(HE_0\omega )^2]_{\text{eff}}=O^{}(HE_0\omega )^2O,$$
(21)
where the operator $`O`$ represents the projection onto the superblock subspace. For a typical many-body Hamiltonian $`H`$ such a calculation is excessively complicated and computationally intensive. Therefore, I calculate an effective representation of $`H`$ only, $`H_{\text{eff}}=O^{}HO`$, and assume that

$$[(HE_0\omega )^2]_{\text{eff}}(H_{\text{eff}}E_0\omega )^2$$
(22)
to calculate $`W_{A,\eta }(\omega ,\psi )`$ in the superblock subspace. This second approximation can cause a violation of the variational bound $`W_{A,\eta }(\omega ,\psi )\pi \eta I_A(\omega +i\eta )`$. Fortunately, the substitution (22) has no significant effect on the minimum of $`W_{A,\eta }(\omega ,\psi )`$ if the state $`(HE_0\omega )|Y_A(\omega +i\eta )|X_A(\omega +i\eta )`$ is accurately represented in the DMRG basis \[i.e., if $`O|X_A(\omega +i\eta )|X_A(\omega +i\eta )`$ for all superblock subspaces\]. Therefore, to use the substitution (22) without loss of accuracy it is necessary and sufficient to include the state $`|X_A(\omega +i\eta )`$ as a target in a DDMRG calculation, even if the real part of the dynamical correlation function is not calculated.

In practice, for sufficiently large $`m`$ I have found that the absolute values of errors in a spectrum $`I_A(\omega +i\eta )`$ decrease systematically with increasing $`m`$. Therefore, it is possible to estimate the accuracy of a DDMRG calculation from the results obtained for different values of $`m`$ as one can do for static properties calculated with DMRG. Moreover, DDMRG results for $`I_A(\omega +i\eta )`$ tend to be smaller than the exact result for almost all frequencies although they can exceed it occasionally.

Obviously, DDMRG is very similar to the correction vector DMRG. till ; pati  The same DMRG basis is built because the same target states are used in both methods. As numerical errors are usually dominated by the DMRG basis truncation, the correction vector parts $`|X_A(\omega +i\eta )`$ and $`|Y_A(\omega +i\eta )`$ are calculated with the same precision $`ϵ`$ in both methods for a given number $`m`$ of density-matrix eigenstates kept per block. Nevertheless, the variational formulation has two significant advantages. First, the errors in the spectrum $`I_A(\omega +i\eta )`$, the excitation energies $`E_nE_0`$, and the spectral weights are of the order of $`ϵ^2`$ instead of $`ϵ`$ in the correction vector method, as explained in Sec. II. If one uses the Lanczos algorithm instead of Eq. (8) in the correction vector DMRG, errors become even larger. Therefore, DDMRG results are more accurate than those obtained with the correction vector DMRG for a given number $`m`$ of density-matrix eigenstates. Second, a DDMRG calculation is essentially an application of the standard DMRG algorithm to the minimization of a different functional. In particular, the numerical accuracy and computational effort are controlled by the sole parameter $`m`$ in an optimized DDMRG calculation as in a ground-state DMRG calculation. The correction vector DMRG till ; pati  and Lanczos DMRG karen ; till  are significantly more complicated than the standard DMRG method. In particular, the numerical accuracy and computational effort depend significantly and sometimes unpredictably on the specific states (Lanczos vectors and correction vectors) included as target in the density-matrix renormalization. Therefore, it is easier to implement and use DDMRG than the correction vector DMRG or the Lanczos DMRG.

## IV Finite-size scaling

DDMRG allows us to calculate spectral functions of a correlated electron (or spin) system on a finite lattice with a broadening given by the parameter $`\eta >0`$. They have the generic form

$$I_{N,\eta }(\omega )=\frac{1}{\pi }\underset{n}{}A_n(N)\frac{\eta }{(\omega _n(N)\omega )^2+\eta ^2},$$
(23)
where $`\omega _n(N)`$ denotes the excitation energy and $`A_n(N)>0`$ the spectral weight of the system eigenstates, and $`N`$ is the number of lattice sites. Such spectra are discrete for $`\eta 0`$ because there is only a finite number of eigenstates in a finite system. In the thermodynamic limit a spectral function

$$I(\omega )=\underset{\eta 0}{lim}\underset{N\mathrm{}}{lim}I_{N,\eta }(\omega )$$
(24)
can contain discrete and continuous parts. (It should be noted that the order of limits in the above formula is important.) To determine the properties of a dynamical spectrum $`I(\omega )`$ in the thermodynamic limit one has to analyze the scaling of the corresponding spectra $`I_{N,\eta }(\omega )`$ as a function of system size. Here I present a finite-size scaling technique for spectral functions calculated with a numerical method such as DDMRG.

Computing both limits in Eq. (24) from numerical results for $`I_{N,\eta }(\omega )`$ requires a lot of accurate data for different values of $`\eta `$ and $`N`$ and can be the source of large extrapolation errors. A much better approach is to use a broadening $`\eta (N)>0`$ which decreases with increasing $`N`$ and vanishes in the thermodynamic limit. The dynamical spectrum is then given by

$`I(\omega )`$ $`=`$ $`\underset{N\mathrm{}}{lim}I_{N,\eta (N)}(\omega )`$
$`=`$ $`\underset{N\mathrm{}}{lim}{\displaystyle \frac{1}{\pi }}{\displaystyle \underset{n}{}}A_n(N){\displaystyle \frac{\eta (N)}{(\omega _n(N)\omega )^2+\eta ^2(N)}}.`$
From the existence of both limits in Eq. (24) it can be demonstrated that there exists a minimal broadening $`\eta _0(N)0`$, which decreases as a function of $`N`$ and converges to zero for $`N\mathrm{}`$, such that the above equation is valid for all functions $`\eta (N)`$ with $`\eta (N)>\eta _0(N)`$ and $`lim_N\mathrm{}\eta (N)=0`$. The function $`\eta _0(N)`$ depends on the frequency $`\omega `$ considered. For a finite lattice with $`N`$ sites, let $`M_{\omega ,ϵ}(N)`$ be the number of excited states contributing to the spectral function in a small interval of width $`ϵ`$ around the frequency $`\omega `$ (i.e., $`|\omega _n(N)\omega |<ϵ/2`$). If $`M_{\omega ,ϵ}(N)`$ remains finite for any $`ϵ>0`$ as $`N\mathrm{}`$, the spectral function $`I(\omega )`$ is discrete at $`\omega `$ and $`\eta _0(N)=0`$. Equivalently, one can take the $`\eta 0`$ limit first in Eq. (24). If $`M_{\omega ,ϵ}(N)`$ diverges for all $`ϵ>0`$ as $`N\mathrm{}`$, the spectrum is dense at $`\omega `$ and a minimal broadening $`\eta _0(N)>0`$ is required for Eq. (IV) to be valid. For instance, $`\eta _0(N)`$ must be larger than the distance $`\delta \omega =\omega _{n+1}(N)\omega _n(N)`$ between two consecutive excited states in the spectrum. Note that while a continuous spectrum is obviously dense, a dense spectrum can be continuous or discrete. For instance, an infinite number of excited states with $`A_n(N)>0`$ can converge to the same energy as $`N\mathrm{}`$. This seems to happen for the optical conductivity associated with an exciton in a open chain. fabian

The function $`\eta _0(N)`$ depends naturally on the specific problem studied \[i.e., the scaling of the energies $`\omega _n(N)`$ and spectral weights $`A_n(N)`$\]. For the optical conductivity of one-dimensional correlated electron systems such as the Peierls-Hubbard model, I have found numerically that a sufficient condition for all frequencies $`\omega `$ in a dense part of the optical spectrum is

$$\eta \frac{c}{N},$$
(26)
where the constant $`c`$ is comparable to the width of the dynamical spectrum $`I(\omega )`$, which is finite in such lattice models. Usually, one wants to keep the broadening $`\eta `$ as small as possible because it reduces the resolution of the spectrum. Therefore, I use

$$\eta (N)=\frac{c}{N}$$
(27)
in Eq. (IV) to analyze the finite-size scaling of spectral functions calculated with DDMRG and to extrapolate the finite-size results to the thermodynamic limit. The condition (26) has a very simple physical interpretation. The spectral function $`I_{N,\eta }(\omega )`$ represents the dynamical response of the system over a time period $`1/\eta `$ after one has started to apply an external force. Typically, in a lattice model the spectral width is proportional to the velocity of the excitations involved in the system response. Thus the condition (26) means that excitations are too slow to travel the full length $`N`$ of the system in the time interval $`1/\eta `$ and do not “sense” that the system is finite.

An additional benefit of a broadening satisfying the condition (26) is that the finite-system spectrum $`I_{N,\eta }(\omega )`$ becomes indistinguishable from the infinite-system spectrum with the same broadening $`\eta `$ for relatively small $`N`$,

$$I_{N,\eta }(\omega )C_\eta [I(\omega )].$$
(28)
Therefore, if one knows an exact or conjectured spectral function $`I(\omega )`$ for an infinite system, its convolution with a Lorentzian of width $`\eta `$ can be compared directly with the numerical results for the finite-system spectrum $`I_{N,\eta }(\omega )`$. This approach has been applied successfully to the optical conductivity of one-dimensional Mott insulators in Refs. eric,  and fabian, , and additional examples are presented in the next section.

In practice, the extrapolation scheme (IV) works well at fixed frequency for the continuous parts and the non-dense discrete parts of a spectrum $`I(\omega )`$ only. To detect singularities in $`I(\omega )`$ and determine their properties, it is generally easier to analyze the scaling of maxima in $`I_{N,\eta }(\omega )`$ or in its derivative as $`N\mathrm{}`$ and $`\eta 0`$. To perform this scaling analysis one can use a size-dependent broadening $`\eta (N)`$ such that Eq. (IV) is valid and $`I_{N,\eta }(\omega )`$ is a good approximation of $`C_\eta [I(\omega )]`$ around the maximum. Then the scaling of a maximum in $`I_{N,\eta }(\omega )`$ for $`\eta (N)0`$ gives the scaling of the corresponding maximum in $`C_\eta [I(\omega )]`$ for $`\eta 0`$. Here I discuss some examples of this technique which are useful for the analysis of the Peierls-Hubbard model optical conductivity presented in the next section. \[One should also note that to detect the presence of a gap between two bands in an infinite-system spectrum $`I(\omega )`$, it is often faster and more reliable to investigate the scaling of minima in $`I_{N,\eta }(\omega )`$ for $`\eta (N)0`$ than to perform extrapolations at fixed frequencies.\]

First, we consider an infinite-system spectrum with a peak in a continuous band,

$$I(\omega )=I_0\delta (\omega \omega _0)+I_{\text{cont}}(\omega )$$
(29)
for $`|\omega \omega _0|<\mathrm{\Lambda }`$, where $`I_0>0`$ and $`I_{\text{cont}}(\omega )`$ is a continuous function. It is easy to show that for $`\eta \mathrm{\Lambda }`$ the maximum of $`C_\eta [I(\omega )]`$ diverges as $`I_0/(\pi \eta )`$ and that the position of the maximum converges to $`\omega _0`$ for $`\eta 0`$. The maximum of the corresponding finite-system spectra $`I_{N,\eta }(\omega )`$ has the same scaling properties for $`\eta (N)0`$. Therefore, it is possible to detect such a $`\delta `$-peak and determine its weight $`I_0`$ in an infinite-system spectrum, even if $`I_0`$ is only a small fraction of the total spectral weight of $`I_{\text{cont}}`$.

Second, we consider an infinite-system spectrum with a power-law divergence at the band edge

$$I(\omega )=I_0\theta (\omega \omega _0)|\omega \omega _0|^\alpha $$
(30)
for $`|\omega \omega _0|<\mathrm{\Lambda }`$, where $`\theta (x)=0`$ for $`x<0`$ and $`\theta (x)=1`$ for $`x>0`$, $`I_0>0`$, and $`0<\alpha <1`$. One can show that for $`\eta \mathrm{\Lambda }`$ the maximum of $`C_\eta [I(\omega )]`$ diverges as $`\eta ^\alpha `$ and that the position of the maximum converges to $`\omega _0`$ from above for $`\eta 0`$. Again, the maximum of the corresponding finite-system spectra $`I_{N,\eta }(\omega )`$ has the same scaling properties for $`\eta (N)0`$. Thus it is possible to detect such a singularity and determine the exponent $`\alpha `$ from the finite-system numerical data.

Third, we consider a continuous infinite-system spectrum with a singularity in its derivative

$$I(\omega )=I_0\theta (\omega \omega _0)|\omega \omega _0|^\alpha $$
(31)
for $`|\omega \omega _0|<\mathrm{\Lambda }`$, where the function $`\theta (x)`$ and the constants $`I_0`$ and $`\alpha `$ are as in the previous example. Here the derivative of $`C_\eta [I(\omega )]`$ has a maximum which diverges as $`\eta ^{\alpha 1}`$ for $`\eta \mathrm{\Lambda }`$ while its position converges to $`\omega _0`$ from above as $`\eta 0`$. The maximum in the derivative of $`I_{N,\eta }(\omega )`$ has the same scaling properties for $`\eta (N)0`$. Therefore, it is possible to determine the exponent $`\alpha `$ from the finite-system numerical data in this case too.

Finally, we consider a special function representing a continuous spectrum above a gap and a truncated divergence close to the band edge

$$I(\omega )=I_0\theta (\omega \omega _0)\frac{2\sqrt{\omega _0|\omega \omega _0|}}{\gamma \omega _0+|\omega \omega _0|}$$
(32)
for $`|\omega \omega _0|<\mathrm{\Lambda }`$, where $`\gamma `$ is a constant such that $`0\gamma <\mathrm{\Lambda }/\omega _0`$. The function $`\theta (x)`$ and the other constants are as in the previous examples. For $`\gamma >0`$ this spectrum vanishes as a square-root at the band edge $`\omega _0`$, goes trough a maximum $`I_0/\sqrt{\gamma }`$ at $`\omega =(1+\gamma )\omega _0`$ then decreases monotonically. For $`\gamma 1`$ the maximum is very sharp and close to the band edge, and $`I(\omega )`$ appears to diverge as $`1/\sqrt{|\omega \omega _0|}`$ at higher energy $`\omega >(1+\gamma )\omega _0`$. The continuous vanishing of $`I(\omega )`$ at the band edge is apparent only in a small frequency range $`\omega _0<\omega <(1+\gamma )\omega _0`$. As $`\gamma 0`$ this maximum becomes a square-root singularity at $`\omega _0`$. A qualitatively similar behavior is often found in the optical conductivity of one-dimensional insulators (see the discussion in the next section). Obviously, the maximum of $`C_\eta [I(\omega )]`$ tends to $`I_0/\sqrt{\gamma }`$ and its position converges to $`\omega =(1+\gamma )\omega _0`$ for $`\eta 0`$. For $`\gamma <1`$, however, the convergence of the maximum becomes apparent only for $`\eta \gamma \omega _0`$. For larger $`\eta `$ the maximum appears to diverge as $`1/\sqrt{\eta }`$. Similarly, the maximum of the derivative diverges as $`1/\sqrt{\eta }`$ for $`\eta 0`$ as discussed the previous example. In the present case, however, this scaling is not observed as soon as $`\eta \mathrm{\Lambda }`$ but only if $`\eta \gamma \omega _0`$. The finite-system spectrum $`I_{N,\eta }(\omega )`$ and its derivative have the same scaling properties for $`\eta (N)0`$. Therefore, with the scaling analysis of finite-system spectra it is possible to distinguish an infinite-system spectrum with a truncated divergence above the band edge from a spectrum with a real divergence at the band edge, provided that one can do calculations with a resolution $`\eta (N)\gamma \omega _0`$.

In summary, the dynamical spectrum of an infinite system can be determined accurately and efficiently from numerical data for finite-system spectra using a size-dependent broadening $`\eta (N)`$. The broadening $`\eta (N)`$ must be larger than a minimal broadening $`\eta _0(N)`$, which depends on the system investigated and can vary with the frequency. Often this broadening conceals the finite-size effects and one can directly compare finite-system spectra to analytical results for infinite systems using a convolution with a Lorentzian distribution, see Eq. (28). If this comparison is not possible or not sufficient, specific points of the spectrum can be extrapolated to the thermodynamic limit using Eq. (IV). Finally, the scaling of maxima in finite-system spectra or their derivatives \[for $`\eta (N)0`$\] allows us to find and analyze singularities in the infinite-system spectrum. For one-dimensional correlated electron systems a sufficient condition for the minimal broadening is given by Eq. (26) and one can use a size-dependent broadening (27).

## V Optical conductivity of the Peierls-Hubbard model

In this section I apply the DDMRG method and the finite-size scaling analysis to the optical conductivity of the one-dimensional Peierls-Hubbard model. dionys  This model is defined by the Hamiltonian

$$H=T+U\underset{l=1}{\overset{N}{}}\left(n_{l,}\frac{1}{2}\right)\left(n_{l,}\frac{1}{2}\right)$$
(33)
with

$$T=\underset{l;\sigma }{}\left(t(1)^l\frac{\mathrm{\Delta }}{2}\right)\left(c_{l,\sigma }^+c_{l+1,\sigma }+c_{l+1,\sigma }^+c_{l,\sigma }\right).$$
(34)
It describes electrons with spin $`\sigma =,`$ which can hop between neighboring sites in a lattice with an even number $`N`$ of sites. In Eq. (34) the index $`l`$ runs from 1 to $`N1`$ for an open chain and from 1 to $`N`$ if periodic boundary conditions are used. Here $`c_{l,\sigma }^+`$, $`c_{l,\sigma }`$ are creation and annihilation operators for electrons with spin $`\sigma `$ at site $`l`$ and $`n_{l,\sigma }=c_{l,\sigma }^+c_{l,\sigma }`$ are the corresponding density operators. The hopping integral $`t>0`$ gives rise to a single-electron band of width $`4t`$. The dimerization parameter $`0|\mathrm{\Delta }|2t`$ determines the strength of the periodic lattice potential generated by the Peierls instability. (For a finite open chain I only use $`\mathrm{\Delta }0`$ to avoid spurious excitations at the chain ends.) The Coulomb repulsion is mimicked by a local Hubbard interaction $`U0`$. The number of electrons equals the number of lattice sites.

The ground state, single-particle charge gap and spin gap of this system can be calculated with great accuracy on lattices with up to $`N10^3`$ sites using DMRG. The single-particle charge gap is given by

$$E_c(N)=E_0(N+1)+E_0(N1)2E_0(N),$$
(35)
where $`E_0(M)`$ denotes the ground-state energy for $`M`$ electrons in a $`N`$-site system. For an even number $`N`$ of sites and electrons the Peierls-Hubbard model ground state is a singlet lieb  and the spin gap is given by

$$E_s(N)=E_0(S_z=\mathrm{})E_0(S_z=0),$$
(36)
where $`S_z`$ is the the $`z`$-component of the total spin and $`E_0(S_z)`$ is the ground-state energy for a fixed value of $`S_z`$.

Spectroscopy with electromagnetic radiation is a common experimental probe of solid-state materials. The linear optical absorption is proportional to the real part $`\sigma _1(\omega )`$ of the optical conductivity. For $`\omega >0`$, $`\sigma _1(\omega )`$ is related to the imaginary part of the current-current correlation function $`G_J(\mathrm{}\omega +i\eta )`$ by

$`\sigma _1(\omega )`$ $`=`$ $`{\displaystyle \frac{\pi }{Na\omega }}\underset{\eta 0}{lim}\text{Im}G_J(\mathrm{}\omega +i\eta )`$
$`=`$ $`{\displaystyle \frac{\pi }{Na\omega }}{\displaystyle \underset{n}{}}|\psi _0|J|\psi _n|^2\delta (\mathrm{}\omega +E_0E_n).`$
Here $`|\psi _0`$ is the ground state of the Hamiltonian $`H`$, $`|\psi _n(n>1)`$ are the other eigenstates of $`H`$, and $`E_0`$, $`E_n`$ are their respective eigenenergies. In this model the current operator is

$$J=\frac{iae}{\mathrm{}}\underset{l;\sigma }{}\left(t(1)^l\frac{\mathrm{\Delta }}{2}\right)\left(c_{l,\sigma }^+c_{l+1,\sigma }c_{l+1,\sigma }^+c_{l,\sigma }\right),$$
(38)
where $`a`$ is the lattice constant, $`e`$ is the charge of an electron, and the index $`l`$ takes the same values depending on the boundary conditions as in Eq. (34). Note that this is the natural definition of the current operator for both types of boundary conditions. The Parzen filter used for open boundary conditions in other works till2 ; brune  is not necessary and thus not used in this work.

In an open chain the optical absorption is also related to the dynamical polarizability $`\alpha (\omega )`$, which is given by the imaginary part of the dipole-dipole correlation function $`G_D(\mathrm{}\omega +i\eta )`$,

$`\alpha (\omega )`$ $`=`$ $`{\displaystyle \frac{\pi }{Na}}\underset{\eta 0}{lim}\text{Im}G_D(\mathrm{}\omega +i\eta )`$
$`=`$ $`{\displaystyle \frac{\pi }{Na}}{\displaystyle \underset{n}{}}|\psi _0|D|\psi _n|^2\delta (\mathrm{}\omega +E_0E_n).`$
For a constant lattice spacing $`a`$ the dipole operator is

$$D=ea\underset{l=1}{\overset{N}{}}l\left(n_l1\right)$$
(40)
with $`n_l=n_{l,}+n_{l,}`$. Using the relation $`J=\frac{i}{\mathrm{}}[D,H]`$ one easily proves that

$$\sigma _1(\omega )=\omega \alpha (\omega )$$
(41)
and

$$\sigma _1(\omega )=\frac{\pi }{Na\mathrm{}}\underset{\eta 0}{lim}\text{Im}\{(\mathrm{}\omega +i\eta )G_D(\mathrm{}\omega +i\eta )\}.$$
(42)
The optical conductivity can be calculated with the DDMRG method described in this paper. For an open chain Eqs. (V),(V), and (42) provide us with three different approaches. First, one can calculate the imaginary part of the current-current correlation function with DDMRG and use Eq. (V) to obtain the convolution of the reduced optical conductivity

$`C_\eta [\omega \sigma _1(\omega )]`$ $`=`$ $`{\displaystyle \frac{\pi }{Na}}\text{Im}G_J(\mathrm{}\omega +i\eta )`$
$`=`$ $`{\displaystyle \frac{1}{Na}}{\displaystyle \underset{n}{}}{\displaystyle \frac{\eta |\psi _0|J|\psi _n|^2}{(\mathrm{}\omega +E_0E_n)^2+\eta ^2}}.`$
This is also the only approach possible with periodic boundary conditions. Second, one can calculate the imaginary part of the dipole-dipole correlation function with DDMRG and use Eq. (V) to obtain the convolution of the dynamical polarizability

$`C_\eta [\alpha (\omega )]`$ $`=`$ $`{\displaystyle \frac{\pi }{Na}}\text{Im}G_D(\mathrm{}\omega +i\eta )`$
$`=`$ $`{\displaystyle \frac{1}{Na}}{\displaystyle \underset{n}{}}{\displaystyle \frac{\eta |\psi _0|D|\psi _n|^2}{(\mathrm{}\omega +E_0E_n)^2+\eta ^2}}.`$
The optical conductivity is then given by the relation (41). Finally, one can calculate the complete dipole-dipole correlation function with DDMRG and use Eq. (42) to obtain the convolution of $`\sigma _1(\omega )`$ directly

$`C_\eta [\sigma _1(\omega )]`$ $`=`$ $`{\displaystyle \frac{\pi }{Na\mathrm{}}}\text{Im}\{(\mathrm{}\omega +i\eta )G_D(\mathrm{}\omega +i\eta )\}`$
$`=`$ $`{\displaystyle \frac{1}{Na\mathrm{}}}{\displaystyle \underset{n}{}}{\displaystyle \frac{\eta |\psi _0|D|\psi _n|^2(E_nE_0)}{(\mathrm{}\omega +E_0E_n)^2+\eta ^2}}.`$
$`C_\eta [\sigma _1(\omega )]`$ can also be formulated in terms of the current matrix elements $`|\psi _0|J|\psi _n|^2`$ \[see Eq. (4) of Ref. fabian, \]. Although the real part of a dynamical correlation function is used in Eq. (V) to calculate the optical conductivity, its relative contribution to $`C_\eta [\sigma _1(\omega )]`$ is of the order $`(\eta /t)^2`$. Therefore, the numerical precision is not significantly reduced by the lower accuracy of DDMRG for the real part of dynamical correlation functions.

Clearly, all three approaches give the same spectrum $`\sigma _1(\omega )`$ for $`\eta 0`$. In DDMRG calculations with $`\eta >0`$, however, they are not equivalent. First, I have found that it is easier to calculate the dipole-dipole correlation function $`G_D(\mathrm{}\omega +i\eta )`$ than the current-current correlation function $`G_J(\mathrm{}\omega +i\eta )`$, except for very strong coupling $`Ut`$. Second, the finite-size scaling \[using a size-dependent broadening (27)\] is different for the three optical spectra $`C_\eta [\omega \sigma _1(\omega )]/\omega `$, $`\omega C_\eta [\alpha (\omega )]`$, and $`C_\eta [\sigma _1(\omega )]`$, especially for very small and very large frequencies $`\omega `$. Usually, $`C_\eta [\sigma _1(\omega )]`$ is the best approximation to $`\sigma _1(\omega )`$ but at low energy ($`\mathrm{}\omega <t`$) it can be more convenient to use $`\omega C_\eta [\alpha (\omega )]`$ while at high energy ($`\mathrm{}\omega t`$) I prefer $`C_\eta [\omega \sigma _1(\omega )]`$.

To calculate the optical spectrum of the Peierls-Hubbard model I have used the third approach, Eq. (V), in most cases. Thus $`C_\eta [\sigma _1(\omega )]`$ is shown in the figures of this paper unless I state explicitly otherwise. Only optical spectra calculated with open boundary conditions are presented. As the size-dependent broadening (27) conceals most of the finite-size effects, spectra calculated with periodic boundary conditions would be almost identical. In all figures showing optical spectra I set $`a=e=\mathrm{}=t=1`$. Thus $`\sigma _1(\omega )`$ is shown in units of $`e^2a/\mathrm{}`$, $`\omega \sigma _1(\omega )`$ in units of $`e^2at/\mathrm{}^2`$, and the frequency $`\omega `$ in units of $`\mathrm{}/t`$.

The sum rules (6) take a simple form for the optical conductivity in the Peierls-Hubbard model with open boundary conditions

$`{\displaystyle \frac{\mathrm{}}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega \omega \sigma _1(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{Na}}\psi _0|J^2|\psi _0,`$ (46a)
$`{\displaystyle \frac{\mathrm{}}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega \sigma _1(\omega )`$ $`=`$ $`{\displaystyle \frac{ae^2}{2N\mathrm{}}}\psi _0|T|\psi _0,`$ (46b)
$`{\displaystyle \frac{\mathrm{}}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑\omega {\displaystyle \frac{\sigma _1(\omega )}{\omega }}`$ $`=`$ $`{\displaystyle \frac{1}{Na}}\psi _0|D^2|\psi _0.`$ (46c)
To prove the second sum rule (46b) one uses the relation $`[D,J]=ia^2e^2T/\mathrm{}`$. With periodic boundary conditions, only the first two sum rules remain valid. In the second sum rule, however, one must take into account the coherent part of the conductivity at $`\omega =0`$ and the proof is more complicated than for an open chain.

The right-hand-side of Eq. (46) can be calculated accurately with the ground-state or Lanczos DMRG method. Note that for $`\eta >0`$ an optical spectrum calculated from Eq. (V), (V), or (V) exactly fulfills the sum rule (46a), (46b), or (46c), respectively. For the DDMRG results presented in this paper, the sum rules are fulfilled within a few percents.

The optical gap $`E_{\text{opt}}(N)`$ is defined as the excitation energy ($`E_nE_0`$) of the lowest eigenstate contributing to the optical conductivity (i.e., $`\psi _n|J|\psi _00`$) in a $`N`$-site system. $`E_{\text{opt}}(N)`$ can be calculated with the DDMRG method for individual excited states described in Sec. III. As the Peierls-Hubbard Hamiltonian (33) has a particle-hole symmetry the optical gap can also be determined using the symmetrized DMRG method. ramasesha  As expected, both approaches give the same results for $`E_{\text{opt}}(N)`$ within numerical errors. In the thermodynamic limit ($`N\mathrm{}`$) I have found that the optical gap $`E_{\text{opt}}`$ is equal to the single-particle charge gap $`E_c`$ (35) for all $`U0`$ and $`2t>\mathrm{\Delta }0`$. \[In the dimer limit ($`\mathrm{\Delta }=2t`$) the Hamiltonian (33) describes independent dimers: the optical weight is concentrated in a single peak corresponding to Frenkel excitons localized on a dimer and the other (delocalized) excitations above and below this peak carry not optical weight, thus $`E_c<E_{\text{opt}}`$.\] In a finite system or with additional electronic interactions fabian  the single-particle charge gap $`E_c`$ can be different from the optical gap $`E_{\text{opt}}`$.

All DMRG methods have a truncation error which is reduced by increasing the number $`m`$ of retained density matrix eigenstates (for more details, see Refs. steve,  and dmrgbook, ). Varying $`m`$ allows one to compute physical quantities for different truncation errors and thus to obtain error estimates on these quantities. For some quantities, especially eigenenergies, it is possible to extrapolate the results to the limit of vanishing truncation error and thus to achieve a greater accuracy. I have systematically used these procedures to estimate the precision of my numerical calculations and adjusted the maximal number $`m`$ of density matrix eigenstates to reach a desired accuracy. This is especially important for DDMRG calculations as truncation errors in dynamical spectra can greatly vary as a function of the frequency $`\omega `$ for fixed $`m`$. In this work the largest number of density matrix eigenstates used is $`m=600`$. For all numerical results presented here DMRG truncation errors are negligible.

In the following three subsections I demonstrate the finite-size scaling technique and the accuracy of DDMRG on three special limits of the Peierls-Hubbard model. Then the optical conductivity of a Mott-Peierls insulator is presented and discussed in the last subsection.

### V.1 Peierls insulator

For $`U=0`$ the Hamiltonian (33) describes a system of independent electrons, which can be solved exactly for any value of $`\mathrm{\Delta }`$, boundary conditions, or lattice size. This provides us with a perfect test case for the DDMRG method. I have checked that DDMRG can reproduce the optical spectrum of this system on lattices with several hundred sites, for any frequency $`\omega `$, and with relative errors as small as $`10^4`$ using only a few hundred density-matrix eigenstates. This demonstrates that one can obtain almost exact results for the optical conductivity of finite one-dimensional systems such as the Peierls-Hubbard model using DDMRG.

In the thermodynamic limit the Hamiltonian (33) describes a Peierls (band) insulating phase for $`\mathrm{\Delta }0`$ and $`U=0`$. The optical gap $`E_{\text{opt}}`$, the charge gap $`E_c`$, and the spin gap $`E_s`$ equal $`2|\mathrm{\Delta }|`$. Optical excitations are made of one hole in the valence band and one electron in the conduction band. The optical conductivity is given by

$$\sigma _1(\omega )=\frac{ae^2(2\mathrm{\Delta })^2(4t)^2}{2\mathrm{}(\mathrm{}\omega )^2\sqrt{[(\mathrm{}\omega )^2(2\mathrm{\Delta })^2][(4t)^2(\mathrm{}\omega )^2]}}$$
(47)
for $`2|\mathrm{\Delta }|<\mathrm{}\omega <4t`$ and is zero elsewhere. florian  This optical spectrum contains a single band of width $`4t2|\mathrm{\Delta }|`$ with square-root divergences at both band edges. These divergences are a typical feature of a one-dimensional band insulator. The convolution of Eq. (47) with a Lorentzian distribution of width $`\eta /t=0.05`$ is shown in Fig. 1 for $`\mathrm{\Delta }=0.6t`$. Both divergences are replaced by maxima at $`\mathrm{}\omega 2\mathrm{\Delta }=1.2t`$ and $`\mathrm{}\omega 4t`$. In Fig. 1 I also show the optical conductivity calculated with DDMRG on a 128-site lattice with the same broadening. We see that the finite-system optical spectrum is indistinguishable from the infinite-system spectrum. The broadening $`\eta /t=0.05`$ satisfies the condition (26) and thus conceals the finite-size effects as discussed in Sec IV. In this case a broadening $`\eta (N)/t=6.4/N`$ is enough because the spectrum band width is smaller than $`4t`$.

With this size-dependent broadening $`\eta (N)`$ one can use Eq. (IV) to extrapolate the finite-size DDMRG results to the thermodynamic limit. For instance, for $`\mathrm{}\omega =2.6t`$ I have obtained $`\sigma _1(\omega )=0.245`$ (in units of $`ae^2/\mathrm{}`$) using data for systems with up to $`N=256`$ sites \[i.e., with a broadening down to $`\eta (N)/t=0.025`$\], in excellent agreement with the exact result 0.243. If we did not know the exact result (47), we could nevertheless determine the existence of square-root divergences at both band edges using a scaling analysis of the maxima in the DDMRG spectra. For instance, the height of the low-energy maximum (close to $`\mathrm{}\omega =1.2t`$) diverges as $`1/\sqrt{\eta }`$ for $`\eta (N)0`$ \[see Fig. 2(a)\]. Moreover, the position of the maximum tends from above to the optical gap $`E_{\text{opt}}=2\mathrm{\Delta }=1.2t`$ for $`N\mathrm{}`$ \[see Fig. 2(b)\]. As explained in Sec. IV these scaling properties correspond to a square-root divergence at the band edge. Figure 2(b) also shows the finite-size optical gaps $`E_{\text{opt}}(N)`$ calculated with the DDMRG method for individual excited states. They tend to the exact result $`E_{\text{opt}}=1.2t`$ for $`N\mathrm{}`$ as expected.

### V.2 Mott-Hubbard insulator

For $`\mathrm{\Delta }=0`$ the Peierls-Hubbard model (33) becomes the one-dimensional Hubbard model at half filling. For $`U>0`$ this model describes a Mott-Hubbard insulator with gapless spin excitations. lieb2  The optical conductivity of this system has recently been determined using DDMRG and analytical methods. eric  Here I only summarize the most important results and give more information about the finite-size scaling analysis carried out in this previous work.

In the half-filled Hubbard model an optical excitation is made of a pair of spinless bosonic excitations carrying opposite charges in the lower (holon) and upper (doublon or anti-holon) Hubbard bands, respectively. As in a Peierls insulator, the optical spectrum consists of a single band but its width is larger, about $`8t`$. A second distinctive feature of this spectrum is a square-root vanishing $`\sigma _1(\omega )\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ at the band threshold $`E_{\text{opt}}`$. There is also a tiny peak in the middle of the band, at least for $`U4t`$.

In Ref. eric,  it is shown that for weak coupling ($`U3t`$) and in the strong-coupling limit ($`U/t\mathrm{}`$) the finite-system optical spectra calculated with DDMRG agree perfectly with the analytical results obtained in the thermodynamic limit using a field-theoretical approach and a strong-coupling analysis, respectively. For instance, Fig. 3 shows the low-energy parts of DDMRG spectra calculated for three different lattice sizes at $`U=3t`$ and the corresponding field-theoretical spectrum for an infinite system. A size-dependent broadening $`\eta (N)/t=12.8/N`$ is used in this case. One clearly sees the convergence of the finite-size spectra toward the field-theoretical result as $`\eta (N)`$ decreases. To make a quantitative comparison one can calculate the convolution of the field-theoretical spectrum with a Lorentzian of width $`\eta `$ satisfying the condition (26) as discussed in Sec. IV. One finds then that finite-size effects are completely concealed by the broadening even for relatively small system sizes. For instance, it is shown in Fig. 3 of Ref. eric,  that the low-energy optical spectrum calculated on a 128-site lattice for $`U=3t`$ is indistinguishable from the field-theoretical spectrum with the same broadening $`\eta /t=0.1`$. In the strong-coupling limit DDMRG and analytical results agree even better and finite-size effects are no longer visible for systems as small as $`N=32`$.

For other coupling strengths ($`4U/t<\mathrm{}`$) it is necessary to analyze the scaling of the finite-system DDMRG spectra to determine the optical conductivity of the Hubbard model in the thermodynamic limit. Using numerical results for lattices with up to $`N=256`$ sites \[i.e., a resolution $`\eta (N)/t`$ down to 0.05\], I have found that for all $`U/t`$ the optical conductivity at the lower band edge has the qualitative behavior described by Eq. (32): $`\sigma _1(\omega )`$ vanishes as $`\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ at the band threshold and there is a maximum in $`\sigma _1(\omega )`$ at a frequency $`\omega =(1+\gamma )E_{\text{opt}}/\mathrm{}`$, where $`\gamma `$ is a small number. Field theory predicts the same behavior with $`\gamma 0.24`$ in the weak-coupling limit while the strong-coupling analysis gives a maximum at $`\mathrm{}\omega =U=E_{\text{opt}}+4t`$, and thus $`\gamma `$ vanishes as $`t/U`$ for $`Ut`$. The distance $`\gamma E_{\text{opt}}`$ between the spectrum threshold and the maximum increases with $`U/t`$. For $`U4t`$ this distance is large enough to determine the finite-size scaling of the lower band edge using systems with up to $`N=256`$ sites. As an example, Fig. 2(a) shows the low-energy maximum in the optical spectrum $`\sigma _1(\omega )`$ calculated with DDMRG for $`U=4t`$ as a function of $`\eta (N)`$. The contrast between the Mott-Hubbard insulator and the Peierls insulator is striking and the maximum in the Mott-Hubbard insulator optical spectrum clearly tends to a constant for $`\eta (N)0`$. \[For $`U=4t`$ the optical gap of the Hubbard model is comparable to that of the Peierls insulator with $`\mathrm{\Delta }=0.6t`$, so that a direct comparison of both systems is relevant.\] In Fig. 2(b) one sees that the finite-size optical gaps calculated with DDMRG converge to the exact result lieb2  $`E_{\text{opt}}=1.287t`$ in the thermodynamic limit, but the maximum tends to a higher energy $`\mathrm{}\omega 1.7t`$. Therefore, one can conclude that there is no divergence at the optical conductivity threshold $`\mathrm{}\omega =E_{\text{opt}}`$. Moreover, it is possible to confirm that $`\sigma _1(\omega )`$ vanishes as $`\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ at the lower band edge using either a similar scaling analysis for the derivative of DDMRG spectra or a direct comparison with the convolution of functions such as Eq. (32) or the field-theoretical optical spectrum. eric

For very weak coupling one would need to calculate $`\sigma _1(\omega )`$ for very large system sizes in order to perform the same scaling analysis. Because $`E_{\text{opt}}`$ vanishes exponentially with $`U/t`$ and the scaling analysis must be performed in the asymptotic regime $`\eta (N)<\gamma E_{\text{opt}}`$, the required system sizes $`N`$ increases exponentially as $`t/E_{\text{opt}}`$ for $`U/t0`$. Fortunately, it is not necessary to carry out this analysis for the Hubbard model because the optical conductivity of the weak-coupling field theory is already in excellent agreement with the optical conductivity of the lattice model for $`U=3t`$.

### V.3 Strong-coupling limit

In this section I discuss the special case of a Mott-Peierls insulator ($`\mathrm{\Delta }0,U>0`$) in the strong-coupling limit $`U/t\mathrm{}`$, for which the shape of the optical spectrum is known analytically. florian2  In this limit there is exactly one electron on each site in the ground state of the Peierls-Hubbard Hamiltonian (33). An optical excitation moves one electron from a site to another and thus creates a double occupation (doublon) and an empty site (holon). Therefore, the optical gap is of the order of $`U`$. These elementary charge excitations are spinless bosons as in the Hubbard model. The properties of the spin degrees of freedom are determined by an effective Heisenberg model with alternating exchange couplings $`J_1(t+\mathrm{\Delta }/2)^2/U`$ and $`J_2(t\mathrm{\Delta }/2)^2/U`$. The spin gap $`E_s`$ vanishes in the limit $`U\mathrm{}`$. However, as there is a gap in the spin excitation spectrum for any finite $`U/t`$ (see also next section), the structure of the spin ground state in the limit $`U/t\mathrm{}`$ is actually similar to that of a gapped state. For instance, the antiferromagnetic spin-spin correlations decreases exponentially with increasing distance. Thus, this strong-coupling limit of the Peierls-Hubbard model is different from the two limiting cases discussed previously and from the general case presented in the next section.

In the thermodynamic limit the optical conductivity can be calculated analytically using some reasonable assumptions. florian2  If $`0<|\mathrm{\Delta }|<2t`$, the spectrum consists of two bands for $`2|\mathrm{\Delta }||\mathrm{}\omega U|2t`$

$$\sigma _1(\omega )=\frac{g_0e^2a}{8\mathrm{}}\frac{\sqrt{[(\mathrm{}\mathrm{\Omega })^2(2\mathrm{\Delta })^2][(4t)^2(\mathrm{}\mathrm{\Omega })^2]}}{\mathrm{}\omega |\mathrm{}\mathrm{\Omega }|},$$
(48)
where $`\mathrm{}\mathrm{\Omega }=\mathrm{}\omega U`$, and a $`\delta `$-peak at $`\mathrm{}\omega =U`$

$$\sigma _1(\omega )=\frac{\pi g_\pi e^2at^2}{\mathrm{}U}\delta (\mathrm{}\omega U)$$
(49)
in the middle of the gap $`4|\mathrm{\Delta }|`$ separating the bands. For $`\mathrm{\Delta }0`$ one recovers the optical spectrum of the Hubbard model in the strong coupling limit, which consists of a single band and a $`\delta `$-peak in the middle of this band. eric ; florian2  The prefactors $`g_0`$ and $`g_\pi `$ are spin form factors given by ground-state spin correlation functions. They are functions of the effective exchange coupling ratio $`J_2/J_1`$ and thus of $`\delta =|\mathrm{\Delta }/2t|`$. Assuming a dimer spin ground state (i.e., $`J_1>0`$ and $`J_2=0`$) one obtains $`g_0=9/4`$ and

$$g_\pi =\frac{1+32\delta +62\delta ^2+32\delta ^3+\delta ^4}{4(1+\delta )^2}.$$
(50)
This result becomes exact in the dimer limit $`|\mathrm{\Delta }|=2t`$, where $`g_\pi =8`$. For $`\mathrm{\Delta }0`$ the dimer spin ground state does not give the correct form factors because it is known exactly that $`g_0+g_\pi =4\mathrm{ln}(2)`$ and it was found numerically that $`g_\pi /g_010^2`$ (see Ref. eric, ).

Figure 4 shows the reduced optical conductivity $`C_\eta [\omega \sigma _1(\omega )]`$ calculated using DDMRG on a 128-site lattice for $`\mathrm{\Delta }=0.6t`$ and $`\eta /t=0.1`$. A logarithmic scale is used to make visible the weak bands on both sides of the strong central peak. In Fig. 4 one can recognize the spectral shape predicted by the strong-coupling analysis. To make a quantitative comparison, however, it is first necessary to determine $`g_0`$ and $`g_\pi `$ using the finite-size scaling analysis of Sec. IV. Here I use a size-dependent broadening $`\eta (N)/t=12.8/N`$ as for the Hubbard model because the spectral width is also of the order of $`8t`$. For $`\mathrm{}\omega =U2t`$ DDMRG results for $`C_\eta [\omega \sigma _1(\omega )]`$ tend to 0.78 (in units of $`e^2at/\mathrm{}^2`$) for $`N\mathrm{}`$. Comparison with Eq. (48) then yields $`g_02.2`$. In Fig. 4 I also show the two bands (48) with this value of $`g_0`$ (without broadening). The agreement with the finite-system DDMRG spectrum is excellent. The small deviations visible close to the band edges are due to the different broadening used for the numerical result ($`\eta /t=0.1`$) and for the analytical result ($`\eta =0`$). They vanish if the same broadening is used in both calculations. Once more this confirms that a broadening satisfying Eq. (26) hides most finite-size effects in this model as already shown by other examples in Sec. V.1 of this paper and in Refs. eric,  and fabian, . In the DDMRG spectra $`C_\eta [\omega \sigma _1(\omega )]`$ the height of the central peak diverges as $`2.42t/\eta `$ (in units of $`e^2at/\mathrm{}^2`$) for increasing $`\eta (N)`$ but its position does not change. This confirms that it corresponds to a $`\delta `$-peak at $`\mathrm{}\omega =U`$ and gives an estimate $`g_\pi 2.42`$. This $`\delta `$-peak broadened with a Lorentzian distribution of width $`\eta /t=0.1`$ is also shown in Fig. 4. One sees that the agreement with the DDMRG result is perfect. A similar finite-size scaling was performed to determine the form factor $`g_\pi `$ in the Hubbard model. eric

One notes the surprisingly good agreement between the form factors determined numerically with DDMRG ($`g_02.2`$ and $`g_\pi 2.42`$) and those obtained using the approximation of a dimer spin ground state ($`g_0=2.25`$ and $`g_\pi 2.52`$). For the value $`\mathrm{\Delta }=0.6t`$ used in this example, the ratio $`J_2/J_10.29`$ of the effective exchange coupling is already quite small and thus the dimer spin ground state is probably a very good approximation of the actual spin ground state.

### V.4 Mott-Peierls insulator

The optical conductivity of the Peierls-Hubbard model is not known for general interaction parameters $`U>0`$ and $`2t>|\mathrm{\Delta }|>0`$. In this regime the system is in a Mott-Peierls insulating phase dionys : both a periodic lattice potential (i.e., the alternating hopping terms) and electronic correlations contribute to the formation of a charge gap $`E_c>0`$ and there is a finite spin gap $`E_s>0`$. Numerical investigations of the charge and spin gaps and of static correlation functions reveal no phase transition at finite $`U`$ and intermediate $`\mathrm{\Delta }`$ (see also Ref. dionys,  and references therein). Thus the entire parameter space ($`0<U/t<\mathrm{}`$, $`0<|\mathrm{\Delta }|<2t`$) belongs to a single Mott-Peierls insulating phase. Figure 5 shows charge and spin gaps as a function of $`U`$ and $`\mathrm{\Delta }`$. These gaps have been calculated on lattices with up to $`N=400`$ sites using DMRG and extrapolated to the thermodynamic limit. The charge gap of the Mott-Peierls insulator is always larger than the gap of the Mott-Hubbard and Peierls insulators in the $`\mathrm{\Delta }=0`$ and $`U=0`$ limits, respectively. The spin gap is always smaller than the charge gap in the Mott-Peierls phase.

In the thermodynamic limit the optical gap $`E_{\text{opt}}`$ is equal to the charge gap. The nature of the optical excitations in the Mott-Peierls insulator is not well understood. Despite the obvious difference between charge and spin excitation energies, $`E_s<E_c`$, it is not even known if there is a spin-charge separation for single-particle excitations. Optical excitations could consist of a pair of fermionic quasi-particles with opposite spins $`\pm \sigma `$ and opposite charges $`\pm e`$ as in a Peierls insulator (Sec. V.1). They could as well be made of two spinless bosonic excitations carrying opposite charges $`\pm e`$ as in the Mott-Hubbard insulator (Sec. V.2) and in the strong-coupling limit (Sec. V.3).

The investigation of spin and charge gaps and static correlation functions clearly shows that the three special cases described in the previous sections are singular limits of the Peierls-Hubbard model. Unsurprisingly, I have found that the optical conductivity in the Mott-Peierls phase is unlike the simple spectrum observed in these limits. \[All optical spectra presented in this section have been calculated using DDMRG and the finite-size scaling analysis has always been performed using a size-dependent broadening $`\eta (N)/t=12.8/N`$.\]

For large but finite $`U`$ the optical spectrum consists of three bands: a narrow band with a strong singularity around $`\mathrm{}\omega =U`$ and one weak band on each side of this central peak. The singularity seems to be made of two very close power-law divergences which merge to form the single isolated $`\delta `$-peak (49) in the $`U/t\mathrm{}`$ limit. The optical spectrum starts as $`\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ at the lower band edge $`E_{\text{opt}}`$ for all $`|\mathrm{\Delta }|<2t`$. Figure 6 shows the reduced optical conductivity $`C_\eta [\omega \sigma _1(\omega )]`$ for $`U=40t`$ and $`\mathrm{\Delta }=0.6t`$. The spectrum looks very similar to the spectrum for $`U/t\mathrm{}`$, which is again shown in this figure. A finite-size-scaling analysis shows however that the strong central peak is not a $`\delta `$-function for $`U=40t`$ but a narrow band with at least one singularity diverging as $`\eta ^{0.8}`$. The spectra in Fig. 6 are made of three bands: the gaps between the bands appear as local minima on both sides of the central peak because of the relatively large broadening used ($`\eta /t=0.1`$). The finite-size scaling analysis confirms the existence of these gaps. For decreasing parameters $`U`$ or $`|\mathrm{\Delta }|`$ first the lower gap, then the upper gap close. Therefore, the number of bands in the optical spectrum of the Mott-Peierls insulator is not constant but depend on the interaction parameters $`U`$ and $`\mathrm{\Delta }`$.

The evolution of the optical conductivity as a function of $`U`$ is very interesting. For decreasing $`U/t`$ one observes that the central peak breaks into two peaks appearing as local maxima in the broadened spectrum of finite-size systems. The first peak (at the lowest energy) takes over most of the optical weight of the central peak. Its weight decreases progressively with decreasing $`U/t`$ but remains strong even for small $`U`$. In Fig. 7, it is clearly visible (at $`\mathrm{}\omega >4t`$) even for $`U=2t`$ (with $`\mathrm{\Delta }=0.6t`$). This peak corresponds to a power-law divergence within a band with an exponent that tends to $`1/2`$ for $`U0`$. The peak position moves to lower energy as $`U`$ decreases and reaches $`\mathrm{}\omega =4t`$ for $`U=0`$. Therefore, the central peak observed at strong coupling $`Ut`$ corresponds to the upper square-root divergence in the Peierls insulator spectrum (47). (In Fig. 7 this divergence is barely visible as a local maximum at $`\mathrm{}\omega 4t`$ because of the relatively large broadening $`\eta /t=0.2`$ used.) The second peak has very little optical weight (it is not visible in Fig. 7) and I have not been able to determine its structure.

For $`Ut`$ the optical weight is distributed symmetrically around the central peak, as seen in Fig. 6. As $`U`$ decreases, there is a progressive transfer of optical weight from high frequency (above the central peak) to low frequency (below the central peak), see Fig. 7. The high-frequency spectrum becomes very weak for small $`U/t`$ but it completely disappears only at $`U=0`$. It is difficult to determine the spectral width for general parameters because the optical conductivity is very weak and vanishes smoothly at high frequency. I estimate that the width of the spectrum lies between $`4t`$ and $`8t`$ for $`U>0`$. The smallest width is reached for large $`\mathrm{\Delta }`$ and small $`U`$, the largest for small $`\mathrm{\Delta }`$ and large $`U`$. The low-frequency spectrum becomes stronger as $`U`$ diminishes. The local maximum below the central peak (seen in Fig. 6) progressively rises, moves closer to the lower band edge, and transforms into a strong narrow peak, visible in the spectra shown in Fig. 7 (at $`\mathrm{}\omega <4t`$). For small enough $`U/t`$ this low-energy peak contains more optical weight than the central peak. For $`U0`$ the low energy peak becomes the square-root divergence of the Peierls insulator spectrum (47) at the band threshold $`E_{\text{opt}}`$.

For $`U>0`$, however, my results suggest that the optical spectrum always vanishes smoothly at the optical gap. I think that the low-energy optical spectrum of the Peierls-Hubbard model at weak coupling has a qualitative behavior similar to that of the Hubbard model: as $`\mathrm{}\omega E_{\text{opt}}`$ decreases, $`\sigma _1(\omega )`$ first appears to diverge as $`(\mathrm{}\omega E_{\text{opt}})^{1/2}`$, then goes through a maximum just above the optical gap $`E_{\text{opt}}`$, and vanishes smoothly for $`\mathrm{}\omega E_{\text{opt}}`$. For large enough $`U`$ it is possible to carry out a finite-size scaling analysis similar to the one performs for the Hubbard model (see Sec. V.2). Thus it is possible to check that the low-energy spectrum maximum is finite and lies at a higher energy than the optical gap, and to show explicitly that $`\sigma _1(\omega )\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ for $`\mathrm{}\omega E_{\text{opt}}0^+`$.

For smaller $`U`$ it becomes increasingly difficult to distinguish a smooth spectrum with a truncated divergence from a true divergence. For instance, Fig. 8 shows the low-frequency optical conductivity $`\omega C_\eta [\sigma _1(\omega )/\omega ]`$ for $`U=2.3t`$ and $`\mathrm{\Delta }=0.15t`$ with a broadening $`\eta /t=0.1`$ ($`N=128`$ sites). For comparison, I also show the spectra in two limits discussed previously, the Mott-Hubbard insulator \[$`\sigma _1(\omega )\sqrt{\mathrm{}\omega E_{\text{opt}}}`$\] and the Peierls insulator \[$`\sigma _1(\omega )1/\sqrt{\mathrm{}\omega E_{\text{opt}}}`$\], with similar optical gaps ($`E_{\text{opt}}/t=0.60.7`$) and the same broadening $`\eta `$. Clearly, the Mott-Peierls insulator spectrum looks like an intermediate case between the spectra observed in both limiting cases. The position of the maximum in the Mott-Peierls insulator spectrum tends to $`0.78t`$ for $`\eta (N)0`$ while the charge gap (and optical gap) equals $`0.704t`$ in the thermodynamic limit. Certainly, there is no divergence in the low-energy optical spectrum. However, the maximum seems to diverge as $`1/\sqrt{\eta }`$ even for the smallest broadening I have used ($`\eta /t=0.05`$). Thus this spectrum seems to be qualitatively similar to a function such as Eq. (32), but the maximum is so close to the optical gap that broadenings $`\eta `$ significantly smaller than $`0.05t`$ (i.e., system sizes much larger than $`N=256`$) would be necessary to reach the asymptotic regime as discussed in Sec. IV. For the same reason it is not possible to determine how the spectrum vanishes for $`\mathrm{}\omega E_{\text{opt}}0^+`$ in such a case

In the Hubbard model it is possible to confirm the absence of a singularity and the square-root vanishing at the band threshold even if the optical gap is as small as $`E_{\text{opt}}=0.4t`$, because we know the optical spectrum of an infinite system for $`E_{\text{opt}}0`$ from field theory. eric  The field theory approach does not only apply to the Hubbard model, but more generally, gives the low-energy optical spectrum of one-dimensional Mott insulators with small Mott gaps. controzzi ; fabian  The different spectral functions depend only on an interaction parameter $`\beta ^21`$. In addition, the optical gap $`E_{\text{opt}}>0`$ and a normalization constant set the frequency scale and the conductivity scale. For $`1/2<\beta ^21`$ these optical spectra described truncated square-root divergence with a square-root vanishing at the band threshold as in Eq. (32). For $`\beta ^2=1`$ (Hubbard model) the spectrum has the shape shown in Fig. 3 with a maximum at $`1.24E_{\text{opt}}`$. As $`\beta ^2`$ decreases the peak becomes sharper and the maximum moves closer to the band edge. For $`\beta ^2=1/2`$ the optical spectrum is similar to that of a Peierls band insulator with a square-root divergence at the band threshold.

Therefore, this field theory fabian ; controzzi  can describe the optical spectrum in both limiting cases (Mott-Hubbard and Peierls insulators) of the Peierls-Hubbard model in the small gap regime, and the field-theoretical spectrum evolves continuously from one limit to the other with $`\beta ^2`$ going from 1 to 1/2. Using $`1\beta ^21/2`$, the optical gap $`E_{\text{opt}}`$, and the normalization constant as fit parameters, I have compared the low-energy optical conductivity calculated with DDMRG for small gaps ($`E_{\text{opt}}<0.71t`$) to field-theoretical spectra with similar broadening as explained in Sec. IV. For instance, I show in Fig. 9 the DDMRG spectrum for the lattice model (33) with $`U=2.3t`$ and $`\mathrm{\Delta }=0.15t`$ and the fitted field-theoretical spectrum with $`\beta ^2=0.58`$. Both spectra agrees up to $`\mathrm{}\omega =1.2t1.7E_{\text{opt}}`$. Generally, I have found that the optical spectrum of a Mott-Peierls insulator can be fitted by a field-theoretical spectrum with $`\beta ^2>1/2`$ over a finite frequency range, from $`\omega =0`$ to a frequency $`\omega `$ which lies between the position of the low-energy maximum and $`2E_{\text{opt}}/\mathrm{}`$. (Naturally, for $`U0`$ the best fit is always obtained with $`\beta ^2=1/2`$.) Therefore, I think that for any $`U>0`$ (and $`0|\mathrm{\Delta }|<2t`$) the optical spectrum vanishes as $`\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ for $`\mathrm{}\omega E_{\text{opt}}0^+`$.

Note that I do not assume that the field-theoretical calculations in Refs. fabian,  and controzzi,  are also valid for the Peierls-Hubbard model with general interaction parameters. Actually, there are visible discrepancies starting at rather low energy between field theory and DDMRG results for the lattice model, as shown in Fig. 9. The agreement between field theory and DDMRG results in the region of the band threshold simply means that the optical spectra in Mott-Peierls insulators and in one-dimensional Mott insulators have similar shapes just above the optical gap.

Finally, it is interesting to examine the evolution of the optical spectrum from weak to strong bond alternation for fixed $`U`$. It has been shown in Ref. eric,  (see also Sec. V.2) that for $`\mathrm{\Delta }=0`$ (Hubbard model) the spectrum consists of a single band with a maximum close to the lower band edge and a tiny peak in the center (at least for $`U4t`$). If $`|\mathrm{\Delta }|`$ increases one observes in Fig. 10(a) that the maximum moves closer to the optical gap and corresponds to a sharper peak. The optical spectrum still starts as $`\sqrt{\mathrm{}\omega E_{\text{opt}}}`$ at the band threshold as discussed above. The central peak, which is too weak to be seen in Fig. 10(a) for $`\mathrm{\Delta }=0`$, becomes rapidly stronger as $`|\mathrm{\Delta }|`$ increases and is clearly visible for $`\mathrm{\Delta }=0.4t`$. As discussed previously this peak becomes a $`\delta `$-function in the strong-coupling limit $`Ut`$ and corresponds to the upper square-root divergence of the Peierls insulator if $`U`$ vanishes. For moderate $`|\mathrm{\Delta }|`$ the ratio between the hopping integrals $`r(\mathrm{\Delta })=(t|\mathrm{\Delta }|/2)/(t+|\mathrm{\Delta }|/2)`$ is not too small and the optical weight is mostly concentrated below the central peak. If this ratio becomes small, however, the central peak becomes the spectrum dominant feature, see Fig. 10(b). The proportion of the optical weight which is in the central peak increases as $`1r^2(\mathrm{\Delta })`$ for $`r(\mathrm{\Delta })0`$. A finite-size scaling analysis confirms however that this peak is not a $`\delta `$-function but is still a power-law divergence within an excitation band. Only in the dimer limit $`|\mathrm{\Delta }|=2t`$ \[$`r(\mathrm{\Delta })=0`$\] the optical spectrum is made of a single $`\delta `$-peak, which corresponds to the excitation of Frenkel excitons localized on a dimer.

In summary, the optical spectrum of a Mott-Peierls insulator consists of one or more bands with a total spectral width ranging from $`4t`$ to $`8t`$. The distinctive features of the spectrum are a square-root vanishing of $`\sigma _1(\omega )`$ at the lower band edge and a peak due to a power-law singularity around the middle of the spectrum. For strong couplings \[$`Ut`$ and $`r(|\mathrm{\Delta }|)1`$\] most of the optical weight is in the central peak, while for weak couplings \[$`U4t`$ and $`r(|\mathrm{\Delta }|)>1/2`$\] it is mostly concentrated in a narrow peak just above the optical gap. In the limit of a vanishing gap ($`U0`$ and $`\mathrm{\Delta }0)`$ this narrow peak becomes a Drude peak at $`\omega =0`$. For intermediate couplings most of the optical weight is distributed over a broad frequency range between the optical gap and the central peak and $`\sigma _1(\omega )`$ goes through a maximum in this range.

The central peak always appears at an energy larger than the bare band width $`4t`$. For parameters which are realistic for conjugated polymers,dionys  most of the optical weight lies below this peak. Therefore, I think that it is not possible to observe such a structure in the optical spectrum of conjugated polymers, because it occurs at a too high energy ($`>`$ 10 eV) and its intensity is too weak.

The optical spectrum of Mott-Peierls insulators is unlike that of Mott-Hubbard and Peierls insulators. However, most of its main features are found in the strong-coupling limit investigated in Ref. florian2,  and discussed in Sec. V.3. This suggests that the optical excitations of a Mott-Peierls insulator could be made of a pair a spinless bosonic excitations with opposite charges as in the strong-coupling limit (and in a Mott-Hubbard insulator). Nevertheless, understanding the nature of the system elementary excitations requires the study of additional dynamical properties, such as the one-particle Green’s functions. The DDMRG method will enable us to carry out this further investigation.

## VI Conclusion

In this paper I have presented a dynamical DMRG method which allows one to calculate the optical conductivity of one-dimensional correlated electron systems on large lattices with great accuracy. The DDMRG approach to the calculation of dynamical properties is essentially an application of the standard DMRG algorithm for ground-state calculations. Therefore, both methods have the same advantages but also the same limitations. In particular, DDMRG will directly benefit from recent and future improvements of DMRG such as the use of additional symmetries.

With DDMRG it is possible to calculate the dynamical response of correlated systems with hundreds of sites and particles. Relative errors of the order of $`10^4`$ can be achieved for the optical spectrum of finite systems with $`N10^2`$ sites. Using a finite-size scaling analysis based on a size-dependent broadening of the discrete finite-system spectra, one can then calculate a dynamical spectrum in the thermodynamic limit with a resolution of the order of 1% of the spectral width and even investigate singularities in a continuum.

The DDMRG approach can be used for various dynamical quantities, such as dynamical spin-spin correlation functions or single-particle Green’s functions. It can also be applied to other lattice quantum many-body models, in higher dimension, including spin or boson degrees of freedom, or long-range interactions. The correction vector DMRG method has been used to calculate non-linear dynamic response functions, such as third-order dynamical polarizabilities. pati  Similarly, the variational principle presented in Sec. II can be generalized to dynamical correlation functions describing these non-linear responses. Thus I believe that it is possible to develop an efficient DDMRG method for calculating these quantities. A limitation of the DDMRG approach is the restriction to zero temperature. It would be desirable to extend the variational principle and the DDMRG approach to finite-temperature dynamical properties.

The computational resources used by the DDMRG method are relatively modest. For instance, all calculations presented in this paper were carried out on workstations with a single 500 MHz Alpha processor and required less than 1 GByte of memory. It would be easy and very efficient to run DDMRG on a parallel computer as calculations for different frequencies are almost independent. This would permit one to investigate much larger or more complicated systems than in this work.

In summary, the DDMRG method and the finite-size scaling technique for dynamical spectra appear extremely accurate and versatile. They provide a powerful approach for investigating the dynamical properties in low-dimensional lattice quantum many-body systems.

Finally, it should be kept in mind that the variational principle for dynamical correlation functions and their related excited states is completely independent from the DMRG method. Therefore, it is possible to combine this principle with other variational methods to calculate dynamical properties. For instance, one could build a trial wavefunction $`|\psi (\{\lambda _i\})=R(\{\lambda _i\})|\psi _0`$, where $`R(\{\lambda _i\})`$ is an operator depending on a few parameters $`\lambda _i`$, such that the calculation of $`W(\{\lambda _i\})=W_{A,\eta }(\omega ,\psi (\{\lambda _i\}))`$ reduces to the evaluation of ground-state correlation functions. Then the minimization of $`W(\{\lambda _i\})`$ with respect to the variational parameters $`\{\lambda _i\}`$ would give a lower bound and an approximate value of the dynamical spectrum $`I_A(\omega +i\eta )`$.

###### Acknowledgements.

I am grateful to F. Essler and F. Gebhard for many stimulating conversations. I also acknowledge useful discussions with T.D. Kühner and S.R. White about the correction vector DMRG method.