text:
# 1 Introduction

## 1 Introduction

Our objective is to extract physics from lattice QCD with possibly minimal amount of computations. Obviously, the required computing power exceeds that of any desktop personal computer currently available in the market. Thus, for one without supercomputer resources, building a computational system  seems to be inevitable if one really wishes to pursue a meaningful number of any physical quantity from lattice QCD. However, the feasibility of such a project depends not only on the funding, but also on the theoretical advancement of the subject, namely, the realization of exact chiral symmetry on the lattice . Now, if we also take into account of the current price/performance of PC hardware components (CPU + RAM + hard disk111The emergence of low-price and high-capacity ($`>100`$ Gbyte) IDE hard disk turns out to be also rather crucial for this project, since the data storage is enormous.), it seems to be the right time to rejuvenate the project  with a new goal - to build a computational system for lattice QCD with exact chiral symmetry. In this paper, we outline the essential features of a Linux PC cluster (64 nodes) which has been built at National Taiwan University. In particular, we discuss how to optimize its hardware and software for lattice QCD with overlap Dirac operator.

First, we start from quenched QCD calculations (i.e., ignoring any internal quark loops by setting $`\text{det}D=1`$). Thus, our first task is to compute quark propagators in the gluon field background, for a sequence of configurations generated stochastically with weight $`\mathrm{exp}(𝒜_g)`$ ($`𝒜_g`$ : pure gluon action). Then the hardronic observables such as meson and baryon correlation functions can be constructed, and from which the hadron masses and decay constants can be extracted. We use the Creutz-Cabbibo-Marinari heat bath algorithm  to generate ensembles of $`SU(3)`$ gauge configurations.

The computation of quark propagators depends on the scheme of lattice fermions, the hard core of lattice QCD. In general, one requires that any quark propagator coupling to physical hadrons must be of the form

$`(D_c+m_q)^1,`$ (1)
where $`m_q`$ is the bare quark mass, and $`D_c`$ is a chirally symmetric and anti-hermitian Dirac operator \[ $`D_c\gamma _5+\gamma _5D_c=0`$ and $`(iD_c)^{}=iD_c`$ \]. Here we assume that $`D_c`$ is doubler-free, has correct continuum behavior, and $`D=D_c(1+raD_c)^1`$ is exponentially local for smooth gauge backgrounds. Note that the way $`m_q`$ coupling to $`D_c`$ is the same as that in the continuum. The chiral symmetry of $`D_c`$ (even at finite lattice spacing) is the crucial feature of any quark coupling to physical hadrons. Otherwise, one could hardly reproduce the low energy strong interaction phenomenology from lattice QCD.

For any massless lattice Dirac operator $`D`$ satisfying the Ginsparg-Wilson relation

$`D\gamma _5+\gamma _5D=2raD\gamma _5D,`$ (2)
it can be written as

$`D=D_c(1+raD_c)^1,`$
and the bare quark mass is naturally added to the $`D_c`$ in the numerator ,

$`D(m_q)=(D_c+m_q)(1+raD_c)^1.`$
Then the quenched quark propagator becomes

$`(D_c+m_q)^1=(1rm_qa)^1[D(m_q)^1ra]`$ (3)
If we fix one of the end points at $`(\stackrel{}{0},0)`$ and use the Hermitcity $`D^{}=\gamma _5D\gamma _5`$, then only 12 (3 colors times 4 Dirac indices) columns of

$`D(m_q)^1=D^{}(m_q)\{D(m_q)D^{}(m_q)\}^1`$ (4)
are needed for computing the time correlation functions of hadrons. Now our problem is how to optimize a PC cluster to compute $`D(m_q)^1`$ for a set of bare quark masses.

The outline of this paper is as follows. In Section 2, we briefly review our scheme of computing quark propagators via the overlap Dirac operator. The details have been given in Ref. . In Section 3, we discuss a simple scheme of memory management for the nested conjugate gradient loops. In Section 4, we discuss how to implement the SSE2 codes for the computationally intense parts of our program. In Section 5, the performance of our system is measured in terms of a number of tests pertaining to the computation of quark propagators. In Section 6, we conclude with some remarks and outlooks.

## 2 Computational Scheme for quark propagators

The massless overlap Dirac operator  reads as

$`D=m_0a^1\left(\text{1I}+\gamma _5{\displaystyle \frac{H_w}{\sqrt{H_w^2}}}\right)`$ (5)
where $`H_w`$ denotes the Hermitian Wilson-Dirac operator with a negative parameter $`m_0`$,

$`H_w=\gamma _5D_w=\gamma _5(m_0+\gamma _\mu t_\mu +W),`$ (6)
$`\gamma _\mu t_\mu `$ the naive fermion operator, and $`W`$ the Wilson term. Then $`D`$ (5) satisfies the Ginsparg-Wilson relation (2) with $`r=1/(2m_0)`$. In this paper, we always fix $`m_0=1.3`$ for our computations. Details of our implementation have been given in Ref. .

Basically, we need to solve the following linear system

$`D(m_q)D^{}(m_q)Y`$ (7)
$`=`$ $`\left\{m_q^2+\left(2m_0^2{\displaystyle \frac{m_q^2}{2}}\right)\left[1+{\displaystyle \frac{(\gamma _5\pm 1)}{2}}H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}\right]\right\}Y=\text{1I}`$
by conjugate gradient (CG). Then the quark propagators can be obtained through (4). With Zolotarev optimal rational approximation  to $`(H_w^2)^{1/2}`$, the multiplication222Note that the Zolotarev optimal rational polynomial in Eq. (8) is in the form $`r^{(n,n)}`$ which is different from $`r^{(n1,n)}`$ used in Ref. . We refer to Ref.  for further discussions.

$`H_w\left({\displaystyle \frac{1}{\sqrt{H_w^2}}}\right)Y,h_w{\displaystyle \frac{H_w}{\lambda _{min}}}`$ (8)
$``$ $`h_w(h_w^2+c_{2n}){\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{b_l}{h_w^2+c_{2l1}}}Y=h_w(h_w^2+c_{2n}){\displaystyle \underset{l=1}{\overset{n}{}}}b_lZ_l`$
can be evaluated by invoking another conjugate gradient process to the linear systems

$`(h_w^2+c_{2l1})Z_l=Y,l=1,\mathrm{},n.`$ (9)
where

$`c_l`$ $`=`$ $`{\displaystyle \frac{\text{sn}^2(\frac{lK^{}}{2n+1};\kappa ^{})}{1\text{sn}^2(\frac{lK^{}}{2n+1};\kappa ^{})}}`$
$`b_l`$ $`=`$ $`d_0{\displaystyle \frac{\underset{i=1}{\overset{n1}{}}(c_{2i}c_{2l1})}{_{i=1,il}^n(c_{2i1}c_{2l1})}}`$
$`d_0`$ $`=`$ $`{\displaystyle \frac{2\lambda }{1+\lambda }}{\displaystyle \underset{l=1}{\overset{n}{}}}{\displaystyle \frac{1+c_{2l1}}{1+c_{2l}}}`$
$`\lambda `$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{2n+1}{}}}{\displaystyle \frac{\mathrm{\Theta }^2(\frac{2lK^{}}{2n+1};\kappa ^{})}{\mathrm{\Theta }^2(\frac{(2l1)K^{}}{2n+1};\kappa ^{})}}.`$
Here $`\mathrm{\Theta }`$ denotes the elliptic theta function, and the Jacobian elliptic function $`\text{sn}(u;\kappa ^{})`$ is defined by the elliptic integral

$`u={\displaystyle _0^{\text{sn}}}{\displaystyle \frac{dt}{\sqrt{(1t^2)(1\kappa ^2t^2)}}},`$
and $`K^{}`$ is the complete elliptic integral of the first kind with modulus $`\kappa ^{}`$,

$`K^{}={\displaystyle _0^1}{\displaystyle \frac{dt}{\sqrt{(1t^2)(1\kappa ^2t^2)}}},`$
where $`\kappa ^{}=\sqrt{11/b}`$, $`b=\lambda _{max}^2/\lambda _{min}^2`$, and $`\lambda _{max}^2`$ and $`\lambda _{min}^2`$ are the maximum and the minimum of the eigenvalues of $`H_w^2`$.

Instead of solving each $`Z_l`$ individually, one can use multi-shift CG algorithm , and obtain all $`Z_l`$ altogether, with only a small fraction of the total time what one had computed each $`Z_l`$ separately. Evidently, one can also apply multi-shift CG algorithm to (7) to obtain several quark propagators with different bare quark masses.

In order to improve the accuracy of the rational approximation as well as to reduce the number of iterations in the inner CG loop, it is crucial to narrow the interval $`[1,b]`$ by projecting out the largest and some low-lying eigenmodes of $`H_w^2`$. We use Arnoldi algorithm  to project these eigenmodes. Denoting these eigenmodes by

$`H_wu_j=\lambda _ju_j,j=1,\mathrm{},k,`$ (10)
then we project the linear systems (9) to the complement of the vector space spanned by these eigenmodes

$`(h_w^2+c_{2l1})\overline{Z_l}=\overline{Y}(1{\displaystyle \underset{j=1}{\overset{k}{}}}u_ju_j^{})Y,l=1,\mathrm{},n.`$ (11)
In the set of projected eigenvalues of $`H_w^2`$, $`\{\lambda _j^2,j=1,\mathrm{},k\}`$, we use $`\lambda _{max}^2`$ and $`\lambda _{min}^2`$ to denote the least upper bound and the greatest lower bound of the eigenvalues of $`\overline{H}_w^2`$, where

$`\overline{H}_w=H_w{\displaystyle \underset{j=1}{\overset{k}{}}}\lambda _ju_ju_j^{}.`$
Then the eigenvalues of

$`h_w^2=\overline{H}_w^2/\lambda _{min}^2`$
fall into the interval $`(1,b)`$, $`b=\lambda _{max}^2/\lambda _{min}^2`$.

Now the matrix-vector multiplication (8) can be expressed in terms of the projected eigenmodes (10) plus the solution obtained from the conjugate gradient loop (11) in the complementary vector space, i.e.,

$`H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}Y{\displaystyle \frac{1}{\lambda _{min}}}H_w(h_w^2+c_{2n}){\displaystyle \underset{l=1}{\overset{n}{}}}b_l\overline{Z}_l+{\displaystyle \underset{j=1}{\overset{k}{}}}{\displaystyle \frac{\lambda _j}{\sqrt{\lambda _j^2}}}u_ju_j^{}YS`$ (12)
Then the breaking of exact chiral symmetry (2) can be measured in terms of

$`\sigma ={\displaystyle \frac{|S^{}SY^{}Y|}{Y^{}Y}}.`$ (13)
In practice, one has no difficulties to attain $`\sigma <10^{12}`$ for most gauge configurations on a finite lattice .

Now the computation of quark propagators involves two nested conjugate gradient loops: the so-called inner CG loop (11), and the outer CG loop (7). The inner CG loop is the price what one pays for preserving the exact chiral symmetry at finite lattice spacing.

## 3 Memory management

In this section we discuss how to configure the hardware and software of a PC cluster such that it can attain the optimal price/performance for the execution of the nested CG loops, (7) and (11).

First, we examine how much memory is required for computing one of the 12 columns of the quark propagators for a set of bare quark masses, since each column can be computed independently. If the required memory can be allocated in a single node, then each node can be assigned to work on one of the 12 columns of the quark propagators. Then the maximum speed of a PC cluster is attained since there is no communication overheads. Nevertheless, the memory (RDRAM) is the most expensive component, thus its amount should be minimized even though the maximum memory at each node can be up to 4 Gbyte. On the other hand, if one distributes the components of the nested CG loops across the nodes and performs parallel computations (with MPI) through a fast network switch, then the memory at each node can be minimal. However, the cost of a fast network switch and its accessories is rather expensive, and also the efficiency of the entire system will be greatly reduced due to the communication overheads. Therefore, to optimize the price/performance of the PC cluster relies on what is the minimal memory required for computing one of the 12 columns of the quark propagators.

Let $`N_s=L^3\times T`$ denote the total number of lattice sites, then each column of $`D^1`$ with double complex (16 bytes) entries takes

$$N_v=N_s\times 12\times 16\text{bytes}.$$
(14)
Using $`N_v`$ or one column as the unit, we list the memory space of all components during the execution of the nested CG loops :

* Gauge links: 3.
* Number of projected low-lying eigenmodes: $`k`$
* Quark propagators \[ i.e., $`Y`$ in (7) \] of $`N_m`$ masses: $`N_m/2`$.
(Note that each $`Y`$ only takes 1/2 column since it is chiral.)
* Conjugate gradient vectors in the CG algorithm: $`N_m/2`$.
* Residual vector for the outer CG loop: $`1/2`$.
* The vector $`\overline{Y}_1`$ (of the smallest bare quark mass) at the interface between the inner and the outer CG loops: 1.
* The inner CG loop: $`2n+3`$ ( where $`n`$ is the degree of Zolotarev rational polynomial), which consists of
(i) $`\{Z_l\}`$ vectors: $`n`$;
(ii) Conjugate gradient vectors $`\{w_l\}`$: $`n`$;
(iii) Residual vector (r): 1;
(iv) $`H_w|w_1`$: 1;
(v) $`H_w^2|w_1`$: 1.

Therefore, the memory space for all components of the nested CG loops is

$`N_{cg}=(N_m+1/2)+(2n+3)+k+3=N_m+2n+k+6.5\text{(columns)}`$ (15)
A schematic diagram of all components of the nested CG loops is sketched in Fig. 1.

Suppose we wish to compute quark propagators on the $`16^3\times 32`$ lattice (at $`\beta =6.0`$), with parameters $`k=16`$, $`n=16`$, and $`N_m=16`$. Then, according to (14) and (15),

$`N_v`$ $``$ $`0.024\text{Gbyte},`$
$`N_{cg}`$ $`=`$ $`70.5\text{columns},`$
the required memory for all components of the nested CG loops is

$`N_{cg}\times N_v70.5\times 0.024=1.7\text{Gbyte}`$
This seems to imply that one should install333At present, most Pentium 4 motherboards designed for housing PC800 RDRAM have 4 memory slots. four stripes of $`512`$ Mbyte modules (i.e. total 2 Gbyte) at each node, if one wishes to let each node compute independently, and to attain the maximum speed of the PC cluster. However, this is a rather expensive solution at this moment, in view of the current price of $`512`$ Mbyte modules. On the other hand, if one distributes the components of the nested CG loops across the nodes and performs parallel computations (with MPI) through a fast network switch, then the price/performance seems to be even worse than the former solution.

Fortunately, we observe that not all column vectors are used simultaneously at any step of the nested CG loops, and also the computationally intense part is at the inner CG loop. Thus we can use the hard disk as the virtual memory for the storage of the intermediate solution vectors and their conjugate gradient vectors ($`Y_\sigma ,P_\sigma ,\sigma =1,\mathrm{}N_m`$) at each iteration of the outer CG loop, while the CPU is working on the inner CG loop. Then the minimal physical memory required at each node can be greatly reduced. Also, the projected eigenmodes are not required to be kept inside the memory, since they are only needed at the start of the inner CG loop to compute $`\overline{Y}_1`$ (for the smallest bare quark mass),

$`\overline{Y}_1(1{\displaystyle \underset{j=1}{\overset{k}{}}}u_ju_j^{})Y_1,`$
and

$`{\displaystyle \underset{j=1}{\overset{k}{}}}{\displaystyle \frac{\lambda _j}{\sqrt{\lambda _j^2}}}u_ju_j^{}Y_1\epsilon _pY_1`$
where $`\epsilon _pY_1`$ is only needed for computing $`S`$ (12) at the completion of the inner CG loop. Thus one has the options to keep the vector $`\epsilon _pY_1`$ inside the memory during the entire inner CG loop or save it to the hard disk and then retrieve it at the completion of the inner CG loop. Further, since $`\overline{Y}_1`$ is only needed at the start of the inner CG loop, so it can share the same memory location with the residual vector $`r`$.

Now it is clear that the minimum memory at each node (without suffering a substantial loss in the performance) is

$`N_{cg}^{min}=(2n+3)+3=2n+6\text{(columns)},`$
which suffices to accommodate the link variables and all relevant vectors for the inner CG loop. After the completion of the inner CG loop and the vector $`S`$ (12) is computed, the memory space of $`2n+3`$ column vectors is released, and the vectors $`\{Y_\sigma \}`$ and $`\{P_\sigma \}`$ of the outer CG loop can be read from the hard disk, which are then updated to new values according to the CG algorithm.

With this simple scheme of memory management, the minimal memory for computing one of the 12 columns of the quark propagators (for a set of bare quark masses) on the $`16^3\times 32`$ lattice with $`n=16`$ (degree of Zolotarev rational polynomial) becomes

$`N_{cg}^{min}\times N_v=38\times 0.024=0.912\text{Gbyte}.`$
Thus the computation can be performed at a single node with one Gbyte of memory, which can be implemented by installing four stripes of 256 Mbyte memory modules, a much more economic solution than using $`4\times 512`$ Mbyte modules. Moreover, the time for disk I/O (at the interface of inner and outer CG loops) only constitutes a few percent of the total time for the execution of the entire nested CG loops (Table 3). This is the optimal memory configuration for a PC cluster to compute quark propagators on the $`16^3\times 32`$ lattice, which of course is not necessarily the optimal one for other lattice sizes. However, our simple scheme of memory management for the nested CG loops should be applicable to any lattice sizes, as well as to other systems.

In passing, we emphasize that the Zolotarev optimal rational approximation to $`(H_w^2)^{1/2}`$ plays a crucial role to minimize the number of vectors required for the inner CG loop. If one had used other rational approximations, then it would require a very large $`n`$ to preserve exact chiral symmetry to a high precision (e.g., $`\sigma <10^{11}`$). In that case, it would be impossible to attain the optimal price/performance as what has been outlined above.

## 4 The SSE2 acceleration

With the optimal memory allocation for each node, we further enhance the performance of our lattice QCD codes (in Fortran) by rewriting its computationally intense parts in the SSE2 assembly codes of Pentium 4. In this section, we briefly review the basic features of the vector unit (SSE2) of Pentium 4, and then describe how to implement SSE2 codes in our lattice QCD program.

### 4.1 The basic features of SSE2

The simplest and the most efficient scheme of parallel computation is Single Instruction Multiple Data (SIMD). It can be implemented inside CPU through a set of long registers. If each register can accommodate several (say, $`s`$) data entries, then any operation (addition, subtraction, multiplication and division) on these registers will act on all data entries in parallel, thus yields the speed-up by a factor of $`s`$ comparing with normal registers. A schematic diagram is shown in Fig. 2.

Even though Intel had implemented the vector unit in their CPUs since Pentium-MMX series, only in the most recent IA-32 Pentium 4 and the advanced IA-64 Itanium, the architecture has been extended to SSE2 (Streamed SIMD Extension 2) to incorporate double precision data entries.

The Pentium 4 processor has eight registers ( %xmm0, %xmm1, …, %xmm7 ) for SIMD operations . Each register is 128 bits wide and can accomodate 4 integers, or 4 single-precision or 2 double-precision floating point numbers. Since we always use double precision floating point numbers in our program, the execution speed of our program can be almost doubled if SSE2 is turned on judiciously in the computationally intensive parts. Note that SSE2 complies with the IEEE 32-bit and 64-bit arithmetic, thus the precision is lower than the extended 80-bit precision of the normal registers in Pentium 4. However, the difference is less than one part in $`10^{15}`$ (double precision), thus is negligible in our computations.

### 4.2 How to implement SSE2 codes in Fortran programs

Since our lattice QCD codes were originally written in Fortran 77, it would be natural if SSE2 codes can be directly embedded in our Fortran program. However, to our knowledge, the Fortran compilers currently available in the market do not support the option of inlining SSE2 codes. Moreover, for optimal performance of SSE2, the data should be aligned to 16-byte memory boundary. This can be easily carried out in C. Therefore our strategy to implement SSE2 codes is rewrite the main program unit in C such that the data arrays are allocated and aligned to 16 bytes memory boundary, then the SSE2 codes are embedded in C subroutines which are then called by original routines in Fortran.

Of course, if one has written lattice QCD codes in C, then the SSE2 codes can be embedded in C routines directly, without dealing with the interface of C and Fortran.

In the following, we illustrate our scheme of implementing SSE2 codes with an example program. The default compilers are $`\mathrm{𝐠𝐜𝐜}`$ and $`\mathrm{𝐠𝟕𝟕}`$ in Linux.

```
        program main
        implicit none
        integer n
        parameter (n=100)
        double precision r(n), v(n)
        call vxzero(n, r, v)
        end

        subroutine vxzero(n, r, v)
        implicit none
        integer n
        double precision c, r(*), v(*)
        ...
        call vadd(n, c, r, v)             !  r = r + c v
        ...                               !  c:    scalar
        end                               !  r, v: vector

```

Here the Fortran main program calls the subroutine vxzero which in turn calls a computationally intensive routine vadd.

First, we rewrite the main program in C, with the data arrays allocated and properly aligned.

```
        #include <malloc.h>
        int main(int argc, char **argv)
        {
            int n=100;
            double *r, *v;
        /* setup the environment for Fortran */
            f_setarg(argc, argv);
            f_setsig();
            f_init();
        /* allocate r & v, and align them to 16-byte boundary */
            r = memalign(16, n*sizeof(double));
            v = memalign(16, n*sizeof(double));
        /* call the Fortran subroutine */
            vxzero_(&n, r, v);
        /* shutdown the I/O channels of Fortran */
            f_exit();
            exit(0);
            return 0;
        }

```

The function call memalign() dynamically allocates 16 bytes aligned pointers r and v. Then the aligned arrays v\[ \] and r\[ \] can be passed to C subroutines for SSE2 operations.

Next we rewrite the computationally intensive routine vadd in C with embedded SSE2 codes.

```
        /* load variable a into %%xmm0 */
        #define sse_load(a)  \
            __asm__ __volatile__ ( "movapd %0, %%xmm0" :: "m" (a))

        /* r = r + %%xmm0 x v          */
        #define sse_add(r, v)               \
            __asm__ __volatile__ (          \
                "movapd %1, %%xmm1 \n\t"    \
                "movapd %2, %%xmm2 \n\t"    \
                "mulpd %%xmm0, %%xmm2 \n\t" \
                "addpd %%xmm1, %%xmm2 \n\t" \
                "movapd %%xmm2, %0"         \
                :                           \
            /* store to address (r), which is indexed as %0 */  \
                "=m" (r)                    \
                :                           \
            /* load from address (r) and (v), which are  \
               indexed as %1 and %2, respectively */     \
                "m" (r), "m" (v))

        #define ALIGN16 __attribute__ ((aligned (16)))

        void vadd_(int *n, double *coeff, double *r, double *v)
        {
            int i, len;
            static double cc[2] ALIGN16;

          /* the array cc is aligned to 16-byte boundary */

            cc[0] = cc[1] = *coeff;
            sse_load(cc[0]);
            len = (*n)/2;
            for (i=0; i<len*2; i+=2) {
                sse_add(r[i], v[i]);
            }
            if (*n % 2 != 0)
                r[len*2] = r[len*2] + cc[0] * v[len*2];
        }

```

Note that we have added the keyword ”\__volatile\__” ( an GNU extension ) in the macro ”\__asm\__”. Its purpose is to ensure that the compiler does not rearrange the order of execution of the codes during compilation. Finally, all object modules are linked by $`\mathrm{𝐠𝐜𝐜}`$ with the option ”-lg2c”.

### 4.3 The implementation of $`H_w`$ times $`|v`$

In our lattice QCD program, most of the execution time is spent in solving quark propagators via the nested CG loops. Thus the execution time is dominated by the operation $`H_w`$ times $`|v`$, which is performed many times ( $`>10^5`$ in most cases ) before the final results of quark propagators can be obtained. Thus it is crucial to optimize this operation with SSE2 codes.

First, we have to set up the correspondence between the data structures used by C and Fortran routines in our program, in particular, for the link variables and the relevant vectors in the nested CG loops.

Suppose we write the arrays of link variables and a column vector $`v`$ in the syntax of Fortran as

$`u(i,j,\mu ,x),`$
$`v(i,k,x),`$
where $`i`$ and $`j`$ are the color indices, $`\mu `$ is the space-time direction, $`k`$ is the spinor index, and $`x`$ is the site index. Now the question is how to access the elements of these arrays in C routines. To resolve this problem, we define some data structures in C as follows.

```
    /* SU(3) matrix, (c01,c02) forms the complex number of u11, and
                     (c03,c04) of u21, etc.                           */
    typedef struct {
        double c01, c02, c03, c04, c05, c06;
        double c07, c08, c09, c10, c11, c12;
        double c13, c14, c15, c16, c17, c18;
    } su3_t;

    /* there are 4 link variables at each site. */
    typedef struct {
        su3_t mu1, mu2, mu3, mu4;
    } ulink_t;

    /* SU(3) vector, (c1,c2) forms the complex number of v1,
                     (c3,c4) of v2, and (c5,c6) of v3. */
    typedef struct {
        double c1, c2, c3, c4, c5, c6;
    } vector_t;

    /* SU(3) Dirac spinor. */
    typedef struct {
        vector_t s1, s2, s3, s4;
    } spinor_t;

```

Then the correspondence can be easily established. For example, the elements $`u(3,2,1,x)`$ and $`v(2,4,x)`$ can be accessed by C routines as ($`u[x]`$.mu1.c11, $`u[x]`$.mu1.c12) and ( $`v[x]`$.s4.c3, $`v[x]`$.s4.c4 ) respectively.

Now we rewrite $`H_w`$ as

$`H_w(x,y)=\gamma _5\left\{(4m_0)\delta _{x,y}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu =1}{\overset{4}{}}}[(1+\gamma _\mu )U_\mu (x)\delta _{x+\mu ,y}(1+\gamma _\mu )U_\mu ^{}(x\mu )\delta _{x\mu ,y}]\right\}`$
Then the multiplication of $`H_w`$ to a column vector $`|v`$ can be optimized by minimizing the number of multiplications involving the link variables. For example, the multiplication in $`(1+\gamma _1)u|v`$ can be written as ( in the spinor space )

$`(\text{1I}+\gamma _1)u|v=\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{cccc}u& 0& 0& u\\ 0& u& u& 0\\ 0& u& u& 0\\ u& 0& 0& u\end{array}\right)\left(\begin{array}{c}v_1\\ v_2\\ v_3\\ v_4\end{array}\right)=\left(\begin{array}{c}u(v_4v_1)\\ u(v_3v_2)\\ r_2\\ r_1\end{array}\right)`$
where all indices are suppressed except the spinor indices. It is clear that the vectors $`v_4v_1`$ and $`v_3v_2`$ should be computed first, before they are multiplied by link variable $`u`$ ( generic symbol for $`U_\mu /2`$ ). For example, the operation $`v_4v_1`$ can be performed by the following macros with SSE2.

```
#define mvpv(v1, v2) \
    __asm__ __volatile__ ( \
        "movapd %0, %%xmm0 \n\t" \
        "movapd %1, %%xmm1 \n\t" \
        "movapd %2, %%xmm2 \n\t" \
        "subpd %3, %%xmm0 \n\t" \
        "subpd %4, %%xmm1 \n\t" \
        "subpd %5, %%xmm2" \
        : : \
        "m" ((v2).c1), \
        "m" ((v2).c3), \
        "m" ((v2).c5), \
        "m" ((v1).c1), \
        "m" ((v1).c3), \
        "m" ((v1).c5))

```

Similarly, we have

$`(\text{1I}\gamma _1)u^{}|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}u^{}(v_4+v_1)\\ u^{}(v_3+v_2)\\ r_2\\ r_1\end{array}\right),`$
$`(\text{1I}+\gamma _2)u|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}u(v_1+iv_4)\\ ir_3\\ u(v_3+iv_2)\\ ir_1\end{array}\right),`$
$`(\text{1I}\gamma _2)u^{}|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}ir_4\\ u^{}(v_2+iv_3)\\ ir_2\\ u^{}(v_4+iv_1)\end{array}\right),`$
$`(\text{1I}+\gamma _3)u|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}u(v_3v_1)\\ u(v_2+v_4)\\ r_1\\ r_2\end{array}\right),`$
$`(\text{1I}\gamma _3)u^{}|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}u^{}(v_1+v_3)\\ u^{}(v_4v_2)\\ r_1\\ r_2\end{array}\right),`$
$`(\text{1I}+\gamma _4)u|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}ir_3\\ ir_4\\ u(v_3+iv_1)\\ u(v_4+iv_2)\end{array}\right),`$
$`(\text{1I}\gamma _4)u^{}|v`$ $`=`$ $`\left(\begin{array}{c}r_1\\ r_2\\ r_3\\ r_4\end{array}\right)=\left(\begin{array}{c}u^{}(v_1+iv_3)\\ u^{}(v_2+iv_4)\\ ir_1\\ ir_2\end{array}\right).`$
So the multiplications involving the link variables can be implemented as

```
    /* for each lattice size */
    for (x=0; x<ldim; x++) {
    /* prefetch for the current multiplication */
        y = iup[x].mu1-1;
        _prefetch_su3(&(u[x].mu1));
        _prefetch_spinor(&(v[y]));
    /* prefetch for the next multiplication */
        z = idn[x].mu1-1;
        _prefetch_su3(&(u[z].mu1));
        _prefetch_spinor(&(v[z]));
    /* r1.s1 = u[x].mu1 * (v[y].s4 - v[y].s1) */
        mvpv(v[y].s1, v[y].s4);
        su3mul(r1.s1, u[x].mu1);
    /* r1.s2 = u[x].mu1 * (v[y].s3 - v[y].s2) */
        mvpv(v[y].s2, v[y].s3);
        su3mul(r1.s2, u[x].mu1);
    /* r1.s3 = -r1.s2 */
        mvset(r1.s3, r1.s2);
    /* r1.s4 = -r1.s1 */
        mvset(r1.s4, r1.s1);

    /* prefetch for the next multiplication */
        y = iup[x].mu2-1;
        _prefetch_su3(&(u[x].mu2));
        _prefetch_spinor(&(v[y]));
    /* r2.s1 = -(u[x].mu1)^{\dagger} * (v[y].s1 + v[y].s4) */
        mvmv(v[z].s1, v[z].s4);
        su3Hmul(r2.s1, u[z].mu1);
    /* r2.s2 = -(u[x].mu1)^{\dagger} * (v[y].s2 + v[y].s3) */
        mvmv(v[z].s2, v[z].s3);
        su3Hmul(r2.s2, u[z].mu1);
    /* r2.s3 = r2.s2 */
        pvset(r2.s3, r2.s2);
    /* r2.s4 = r2.s1 */
        pvset(r2.s4, r2.s1);
        ...

```

where r1, r2, …, and v\[ \] are declared as the type spinor\_t, and u\[ \] is declared as the type ulink\_t. Note that prefetching has been inserted in order to attain the optimal performance. Finally, we have 8 vector segments r1, …, r8, and a diagonal term. They are summed over to give the final result of $`v[y]`$,

$`v[y].s1`$ $`=`$ $`r1.s1+r2.s1+r3.s1+r4.s1+r5.s1+r6.s1+r7.s1+r8.s1`$
$`+(4m_0)v[x].s1,`$
$`v[y].s2`$ $`=`$ $`r1.s2+r2.s2+r3.s2+r4.s2+r5.s2+r6.s2+r7.s2+r8.s2`$
$`+(4m_0)v[x].s2,`$
$`v[y].s3`$ $`=`$ $`(r1.s3+r2.s3+r3.s3+r4.s3+r5.s3+r6.s3+r7.s3+r8.s3`$
$`+(4m_0)v[x].s3),`$
$`v[y].s4`$ $`=`$ $`(r1.s4+r2.s4+r3.s4+r4.s4+r5.s4+r6.s4+r7.s4+r8.s4`$
$`+(4m_0)v[x].s4),`$
Next we come to the question how to implement SSE2 codes for a $`SU(3)`$ matrix times a vector, the most crucial part in $`H_w`$ times $`|v`$. This problem has been solved by Lüscher , and his SSE2 codes is available in the public domain . We found that Lüscher’s code is quite efficient, and have adopted it in our program. For completeness, we briefly outline Lüscher’s algorithm as follows.

Consider

$`\left(\begin{array}{ccc}u_{11}& u_{12}& u_{13}\\ u_{21}& u_{22}& u_{23}\\ u_{31}& u_{32}& u_{33}\end{array}\right)\times \left(\begin{array}{c}y_1\\ y_2\\ y_3\end{array}\right)=\left(\begin{array}{c}r_1\\ r_2\\ r_3\end{array}\right)`$
.
First, the elements ($`y_1,y_2,y_3`$) of the vector $`|y`$ are copied to the registers %xmm0, %xmm1, and %xmm2, respectively. Then the real part of the SU(3) matrix $`\{u_{mn}\}`$ is read sequentially, and is multiplied to $`|y`$ at %xmm0, %xmm1, and %xmm2, and the result is stored at %xmm3, %xmm4, and %xmm5,

$`\%\mathrm{xmm0}`$ $`=`$ $`(\mathrm{Re}(y_1),\mathrm{Im}(y_1)),`$
$`\%\mathrm{xmm1}`$ $`=`$ $`(\mathrm{Re}(y_2),\mathrm{Im}(y_2)),`$
$`\%\mathrm{xmm2}`$ $`=`$ $`(\mathrm{Re}(y_3),\mathrm{Im}(y_3)),`$
$`\%\mathrm{xmm3}`$ $`=`$ $`(t_1,t_2),`$
$`\%\mathrm{xmm4}`$ $`=`$ $`(t_3,t_4),`$
$`\%\mathrm{xmm5}`$ $`=`$ $`(t_5,t_6),`$
where

$`t_1`$ $`=`$ $`\mathrm{Re}(u_{11})\mathrm{Re}(y_1)+\mathrm{Re}(u_{12})\mathrm{Re}(y_2)+\mathrm{Re}(u_{13})\mathrm{Re}(y_3),`$
$`t_2`$ $`=`$ $`\mathrm{Re}(u_{11})\mathrm{Im}(y_1)+\mathrm{Re}(u_{12})\mathrm{Im}(y_2)+\mathrm{Re}(u_{13})\mathrm{Im}(y_3),`$
$`t_3`$ $`=`$ $`\mathrm{Re}(u_{21})\mathrm{Re}(y_1)+\mathrm{Re}(u_{22})\mathrm{Re}(y_2)+\mathrm{Re}(u_{23})\mathrm{Re}(y_3),`$
$`t_4`$ $`=`$ $`\mathrm{Re}(u_{21})\mathrm{Im}(y_1)+\mathrm{Re}(u_{22})\mathrm{Im}(y_2)+\mathrm{Re}(u_{23})\mathrm{Im}(y_3),`$
$`t_5`$ $`=`$ $`\mathrm{Re}(u_{31})\mathrm{Re}(y_1)+\mathrm{Re}(u_{32})\mathrm{Re}(y_2)+\mathrm{Re}(u_{33})\mathrm{Re}(y_3),`$
$`t_6`$ $`=`$ $`\mathrm{Re}(u_{31})\mathrm{Im}(y_1)+\mathrm{Re}(u_{32})\mathrm{Im}(y_2)+\mathrm{Re}(u_{33})\mathrm{Im}(y_3).`$
Next, multiply the vector $`y`$ by $`i=(0,1)`$, i.e.,

$`\begin{array}{ccccc}\%\mathrm{xmm0}\hfill & & (\mathrm{Im}(y_1),\mathrm{Re}(y_1))\hfill & & (\mathrm{Im}(y_1),\mathrm{Re}(y_1)),\hfill \\ \%\mathrm{xmm1}\hfill & & (\mathrm{Im}(y_2),\mathrm{Re}(y_2))\hfill & & (\mathrm{Im}(y_2),\mathrm{Re}(y_2)),\hfill \\ \%\mathrm{xmm2}\hfill & & (\mathrm{Im}(y_3),\mathrm{Re}(y_3))\hfill & & (\mathrm{Im}(y_3),\mathrm{Re}(y_3)),\hfill \end{array}`$
which is implemented by the following SSE2 code

```
    static int sn3[4] ALIGN16 = {0x0,0x80000000,0x0,0x0};
    #define su3mul(r, u) \
        ...                     \
        "xorpd %9, %%xmm0 \n\t" \
        "xorpd %9, %%xmm1 \n\t" \
        "xorpd %9, %%xmm2 \n\t" \
        ...                     \
        ::                      \
        ...                     \
        "m" (sn3[0]));

```

Then the imaginary part of $`\{u_{mn}\}`$ is read and multiplied to $`iy`$, and the final result is

$`\%\mathrm{xmm3}`$ $`=`$ $`(t_1+s_1,t_2+s_2),`$
$`\%\mathrm{xmm4}`$ $`=`$ $`(t_3+s_3,t_4+s_4),`$
$`\%\mathrm{xmm5}`$ $`=`$ $`(t_5+s_5,t_6+s_6),`$
where

$`s_1`$ $`=`$ $`\mathrm{Im}(u_{11})\mathrm{Im}(y_1)\mathrm{Im}(u_{12})\mathrm{Im}(y_2)\mathrm{Im}(u_{13})\mathrm{Im}(y_3),`$
$`s_2`$ $`=`$ $`+\mathrm{Im}(u_{11})\mathrm{Re}(y_1)+\mathrm{Im}(u_{12})\mathrm{Re}(y_2)+\mathrm{Im}(u_{13})\mathrm{Re}(y_3),`$
$`s_3`$ $`=`$ $`\mathrm{Im}(u_{21})\mathrm{Im}(y_1)\mathrm{Im}(u_{22})\mathrm{Im}(y_2)\mathrm{Im}(u_{23})\mathrm{Im}(y_3),`$
$`s_4`$ $`=`$ $`+\mathrm{Im}(u_{21})\mathrm{Re}(y_1)+\mathrm{Im}(u_{22})\mathrm{Re}(y_2)+\mathrm{Im}(u_{23})\mathrm{Re}(y_3),`$
$`s_5`$ $`=`$ $`\mathrm{Im}(u_{31})\mathrm{Im}(y_1)\mathrm{Im}(u_{32})\mathrm{Im}(y_2)\mathrm{Im}(u_{33})\mathrm{Im}(y_3),`$
$`s_6`$ $`=`$ $`+\mathrm{Im}(u_{31})\mathrm{Re}(y_1)+\mathrm{Im}(u_{32})\mathrm{Re}(y_2)+\mathrm{Im}(u_{33})\mathrm{Re}(y_3).`$
## 5 Performance of the system

In this section, we measure the performance of our system by a number of tests pertaining to the computation of quark propagators.

In Table 1, we list the execution time (in unit of second) for $`H_w`$ multiplying a column vector $`Y`$, for both cases with SSE2 turned on and off, and for several lattice sizes. The data shows that turning on SSE2 can speed up our program by a factor $`1.8`$.

In Table 2, we list the execution time (in unit of second) for projecting 20 low-lying eigenmodes of $`H_w^2`$ using ARPACK, versus the number of Arnoldi vectors. It is clear that there exists an optimal number of Arnoldi vectors for a projection, which of course depends on the gauge configuration. In Table 2, the optimal number is $`100`$, which amounts to $`240`$ Mbyte for the $`8^3\times 24`$ lattice. However, for larger lattices such as $`16^3\times 32`$, the optimal number may require more than one gigabyte of memory. In this case, the projection of eigenmodes is carried out at some nodes with 2 gigabyte of memory.

In Table 3, we measure the time used by disk I/O in our simple scheme of memory management for the nested CG loops, versus the number of bare quark masses. The test is performed at a Pentium 4 (2 GHz) node, for the $`16^3\times 32`$ lattice, and with 16 projected eigenmodes. The disk I/O time is the difference of the total execution time between two cases of turning on and off of the memory management. It is remarkable that the percentage of disk I/O time is only $`3\%`$ of the total execution time even for 16 bare quark masses, and with 16 projected eigenmodes. Evidently, for the $`16^3\times 32`$ lattice, our simple scheme of memory management is more efficient and less expensive than any other options, e.g., parallel computing (with MPI) through a fast network switch.

In Table 4, we list the execution time (in unit of second) for a Pentium 4 (2 GHz) node to compute 12 columns of quark propagators in a topologically nontrivial gauge background at $`\beta =5.8`$ on the $`8^3\times 24`$ lattice, versus the number of projected low-lying eigenmodes. Other parameters for the test are: the degree of Zolotarev rational polynomial is $`n=16`$; the number of bare quark masses is $`N_m=12`$; each projected eigenmode satisfies $`(H_w^2\lambda ^2)|x<10^{13}`$, and the stopping criterion for inner and outer CG loops is $`ϵ=10^{11}`$. The execution time is decomposed into three parts : (i) the projections of high and low-lying eigenmodes444Note that the projection time listed in the 2nd column of Table 4 includes 167 seconds for projecting 4 highest eigenmodes of $`H_w^2`$.; (ii) computing 12 columns of $`(DD^{})^1`$ via the nested CG loops; and (iii) computing $`D^{}`$ and multiplying it to $`(DD^{})^1`$. The total time is listed in the last column of the table. For completeness, we also list $`\lambda _{max}`$ and $`\lambda _{min}`$ of $`|\overline{H_w}|`$ (after the projections), the total numbers of iterations of the outer CG loop and average iterations of the inner CG loop, as well as the precision of exact chiral symmetry in terms of $`\sigma `$ (13). Evidently, the time for projecting out the high and low-lying eigenmodes is only a very small fraction of the total execution time for computing 12 columns of quark propagators. However, the projections have very significant impacts on the total execution time since it yields the speed-up by a factor of 2.44, as comparing the first row (no projections) with the last row (projections of 40 low-lying eigenmodes). Moreover, with projections, the exact chiral symmetry can be easily preserved to a very high precision ($`\sigma <10^{13}`$). This suggests that one should project as many low-lying eigenmodes as possible, before executing the nested CG loops. In general, we suspect that the optimal number of projections depends on the projection algorithm, the amount of memory of the system, as well as the gauge configuration.

In Table 5, we measure the precision of exact chiral symmetry $`\sigma `$ (13) versus the degree ($`n`$) of Zolotarev optimal rational polynomial. The values of $`\sigma `$ listed in the second column of Table 5 are the maxima in the nested CG loops. The execution time and the iterations of the nested CG loops are also listed. Evidently, the precision of exact chiral symmetry $`\sigma `$ is quite different from the stopping criterion $`ϵ=10^{11}`$ for inner and outer CG loops, since $`\sigma `$ can be much bigger or smaller than $`ϵ`$, as shown in Table 5, as well as in Tables 4 and 6. It is clear that the necessary condition for preserving exact chiral symmetry to a very high precision is to use a higher degree ($`n`$) Zolotarev rational polynomial for $`(H_w^2)^{1/2}`$. In Ref. , tables are provided for looking up which degree $`n`$ is required to attain one’s desired accuracy in preserving the exact chiral symmetry on the lattice, versus the parameter $`b=\lambda _{max}^2/\lambda _{min}^2`$ of a given gauge configuration.

In Table 6, we list the execution time (in unit of second) of a Pentium 4 (2 GHz) node to compute 12 columns of quark propagators, versus the size of the lattice. The parameters for the test are: the degree of Zolotarev rational polynomial is $`n=16`$, the number of bare quark masses is $`N_m=16`$ and $`ma0.02`$, each projected eigenmode satisfies $`(H_w^2\lambda ^2)|x<10^{13}`$, and the stopping criterion for the inner and the outer CG loops is $`ϵ=10^{11}`$. From the last entry of the last row, we can estimate that a Pentium 4 (2 GHz) node takes about 24 days to complete 12 columns of quark propagators (for 16 bare quark masses) for one gauge configuration at $`\beta =6.0`$ on the $`16^3\times 32`$ lattice. In other words, if we have 12 nodes and let each one of them work on one column of $`D^1`$, then we can complete the quark propagators for one gauge configuration in two days. Since our system consists of 64 nodes, so we can compute quark propagators at a rate more than five gauge configurations per two days.

## 6 Conclusions

In this paper, we outline the essential features of a Linux PC cluster (64 nodes) which has been built at National Taiwan University, and discuss how to optimize its hardware and software for lattice QCD with exact chiral symmetry. At present, all nodes are working around the clock on lattice QCD computations.

With Zolotarev optimal rational approximation to $`(H_w^2)^{1/2}`$, projections of high and low-lying eigenmodes of $`H_w^2`$, the multi-mass CG algorithm, the SSE2 acceleration, and our simple scheme of memory management, we are able to compute quark propagators of $`16`$ bare quark masses on the $`16^3\times 32`$ lattice, with the precision of quark propagators up to $`10^{11}`$ and the precision of exact chiral symmetry up to $`10^{12}`$, at the rate of 2.5 gauge configuration ($`\beta =6.0`$) per day, with our present system of 64 nodes. This demonstrates that an optimized Linux PC cluster can be a viable computational system to extract physical quantities from lattice QCD with exact chiral symmetry .

The speed of our system is higher than 70 Gflops, and the total cost of the hardware is less than US$60000. This amounts to price/performance ratio better than $1.0/Mflops for 64-bit (double precision) computations. The basic idea of optimization is to let each node work independently on one of the 12 columns of the quark propagators (for a set of bare quark masses), and also use the hard disk as the virtual memory for the vectors in the outer CG loop, while the CPU is working on the inner CG loop. Our simple scheme of memory management for the nested CG loops may also be useful to other systems.

In future, we will add more nodes to our system, and will also work on larger lattices, say $`24^3\times 48`$. Then one Gbyte memory at each node is not sufficient to accommodate all relevant vectors in the inner CG loop, even for 12 Zolotarev terms. However, there are several ways to circumvent this problem. First, our memory management scheme is quite versatile, which is more than just for swapping the vectors at the interface of inner and outer CG loops. In fact, it can handle any number of Zolotarev terms for any lattice size, and can automatically minimize disk I/O at any step of the nested CG loops, according to the amount of physical memory of a node. As long as the percentage of the disk I/O time is less than 30%, it is still a better option than distributing the nested CG loops across the nodes and performing parallel computations (with MPI) through a fast network switch, since the communication overheads is expected to be more than 30% of the total time, especially for a system of 100 nodes or more. Secondly, we can increase the amount of memory at each node, which depends on the specification of the motherboard as well as the price and the capacity of the memory modules. Finally, we can also exploit algorithms  which only use five vectors rather than $`2n+3`$ vectors for the inner CG loop, or the Lanczos algorithm as described in Ref. . Now it is clear that a Linux PC cluster is a viable platform to tackle lattice QCD with exact chiral symmetry even for a large lattice ( e.g., $`32^3\times 64`$), though more studies are needed before one reaches an optimal design for dynamical quarks.

Note added in proof:

Recently, it has been shown  that the speed of Neuberger’s double pass algorithm  for computing the matrix-vector product $`R^{(n1,n)}(H_w^2)Y`$ is almost independent of the degree $`n`$ of the rational polynomial, and it is faster than the single pass for $`n>13`$ (for Pentium 4 with SSE2). Thus the single pass has been replaced with the double pass algorithm in our Linux PC cluster.

Acknowledgement

This work was supported in part by the National Science Council, ROC, under the grant number NSC90-2112-M002-021.


format: text

text:
# Optical excitations in a one-dimensional Mott insulator

## I Introduction

In recent years, various quasi-one-dimensional materials, such as conjugated polymers, kiess ; sari  organic charge-transfer salts, farges ; bourbonnais  Cu oxides, hasan  and Ni halides, kishida  have been extensively studied because of their unusual optical properties and their potential application in modern optical technology. In first approximation these materials can be described as one-dimensional strongly correlated electron systems with half-filled bands. The electron-electron interaction drives such a system into a Mott insulating ground state mott ; florianbook  and dominates low-energy excitations. Therefore, the optical properties of one-dimensional Mott insulators are currently a topic of great interest.

Despite many theoretical studies, our knowledge of these systems is still fragmentary because of the difficulties associated with the investigation of strongly correlated systems. For many years, numerical exact diagonalization of small systems was the only method providing reliable information on excited states in correlated electron systems. guo ; diago  Recently, however, the linear optical conductivity and exciton properties of one-dimensional Mott insulators have been calculated analytically for an infinite system in the limit of a small Mott gap JGE ; controzzi ; EGJ  and of a large Mott gap. EGJ ; florian ; gallagher ; stephan ; mizuno ; barford  Moreover, recent developments of the density-matrix renormalization group (DMRG) method steve ; dmrgbook  allows one to calculate excited states and dynamical response functions numerically in large systems and with an accuracy comparable to exact diagonalizations. JGE ; jeckel

A paradigm of a one-dimensional Mott insulator is the extended Hubbard model with hopping integral $`t`$, on-site repulsion $`U`$, and nearest-neighbor repulsion $`V`$ at half filling. Although this model has been widely studied, its properties are still poorly understood in the thermodynamic limit. It is known hirsch  that the system is a Mott insulator in a large region of the parameter space $`(U,V)`$ which is physically relevant for the Coulomb repulsion between electrons $`(U>V0)`$. The precise ground-state phase diagram has only recently be determined using DMRG. jeckel2  In the special case of the Hubbard model hubbard  ($`V=0,U>0`$) the ground state is known exactly to be a Mott insulator lieb  and the optical conductivity has been calculated. JGE  In principle, the optical properties of the Mott insulating phase are also known for $`V>0`$ at weak coupling controzzi ; EGJ  and at strong coupling, EGJ ; florian ; gallagher ; stephan ; mizuno  but the range of validity of these results is not clear a priori. Results for intermediate coupling and close to the Mott phase boundary $`U2V`$ are scarce and the absence of finite-size-effect analysis often hinders their interpretation. stephan ; shuai ; kancharla ; tomita

In this paper I present an accurate and comprehensive investigation of the linear optical conductivity $`\sigma _1(\omega )`$ and the excited states contributing to $`\sigma _1(\omega )`$ in the Mott insulating phase of the one-dimensional extended Hubbard model at half filling. An efficient symmetrized DMRG ramasesha  and the recently developed dynamical DMRG (DDMRG) JGE ; jeckel  are used to calculate the optically excited states and the linear optical conductivity $`\sigma _1(\omega )`$ on large lattices. The comparison of numerical results for $`\sigma _1(\omega )`$ with field-theoretical controzzi ; EGJ  and strong-coupling florian  predictions confirms both the great accuracy of DDMRG and the wide validity range of both analytical methods. I have found that four types of excitations determine the optical properties in the Mott insulating phase: pairs of unbound charge excitations, excitons, excitonic strings, and charge-density-wave (CDW) droplets. I will show that each type of excitations dominate the low-energy spectrum in a particular region of the parameter space $`(U,V)`$ and exhibits a distinct optical spectrum $`\sigma _1(\omega )`$.

The model, the linear optical conductivity $`\sigma _1(\omega )`$, and the relevant symmetries are introduced in detail in the next section. In Sec. III the numerical methods are briefly presented, then the estimation of DMRG truncation errors and finite-size effects are discussed. In Sec. IV I describe the four different types of excitations which contribute to the linear optical conductivity and the corresponding optical spectra in the various interaction regimes from the limit of a large Mott gap to the limit of a small Mott gap. The final section contains the conclusion.

## II Model

The one-dimensional extended Hubbard model is defined by the Hamiltonian

$`\widehat{H}`$ $`=`$ $`t{\displaystyle \underset{l;\sigma }{}}\left(\widehat{c}_{l,\sigma }^+\widehat{c}_{l+1,\sigma }+\widehat{c}_{l+1,\sigma }^+\widehat{c}_{l,\sigma }\right)`$ (1)
$`+U{\displaystyle \underset{l}{}}\left(\widehat{n}_{l,}{\displaystyle \frac{1}{2}}\right)\left(\widehat{n}_{l,}{\displaystyle \frac{1}{2}}\right)`$
$`+V{\displaystyle \underset{l}{}}(\widehat{n}_l1)(\widehat{n}_{l+1}1).`$
It describes electrons with spin $`\sigma =,`$ which can hop between neighboring sites. Here $`\widehat{c}_{l,\sigma }^+`$, $`\widehat{c}_{l,\sigma }`$ are creation and annihilation operators for electrons with spin $`\sigma `$ at site $`l`$, $`\widehat{n}_{l,\sigma }=\widehat{c}_{l,\sigma }^+\widehat{c}_{l,\sigma }`$ are the corresponding density operators, and $`\widehat{n}_l=\widehat{n}_{l,}+\widehat{n}_{l,}`$. The hopping integral $`t>0`$ gives rise to a a single-electron band of width $`4t`$. The Coulomb repulsion is mimicked by a local Hubbard interaction $`U`$, and a nearest-neighbor interaction $`V`$. The physically relevant parameter regime is $`U>V0`$. The number of electrons equals the number of lattice sites $`N`$ (half-filled band). This system is in a Mott insulating phase for $`V<V_c(U)U/2`$ (Ref.hirsch, ). Precise values of the Mott phase boundary $`V_c(U)`$ are given in Ref. jeckel2, .

Note that the chemical potential is chosen in such a way that the Hamiltonian (1) explicitly exhibits a particle-hole symmetry. This Hamiltonian has two other discrete symmetries which are useful for optical excitation calculations: a spin-flip symmetry and a spatial reflection symmetry (through the lattice center). Therefore, each eigenstate has a well-defined parity under charge conjugation ($`P_c=\pm 1`$) and spin flip ($`P_s=\pm 1`$), and belongs to one of the two irreducible representations, $`A_g`$ or $`B_u`$, of a one-dimensional lattice reflection symmetry group.

Spectroscopy with electromagnetic radiation is a common experimental probe of solid-state materials. kuzmany  The linear (one-photon) 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 by

$$\sigma _1(\omega )=\frac{\text{Im}\{\chi (\omega )\}}{\omega }.$$
(2)
For $`\omega 0`$ the current-current correlation function is given by

$`\chi (\omega >0)`$ $`=`$ $`{\displaystyle \frac{1}{Na}}0\left|\widehat{J}{\displaystyle \frac{1}{E_0\widehat{H}+\mathrm{}\omega +i\eta }}\widehat{J}\right|0`$ (3)
$`=`$ $`{\displaystyle \frac{1}{Na}}{\displaystyle \underset{n}{}}{\displaystyle \frac{|0|\widehat{J}|n|^2}{\mathrm{}\omega (E_nE_0)+i\eta }},`$
where $`a`$ is the lattice spacing. Here, $`|0`$ is the ground state of the Hamiltonian $`\widehat{H}`$, $`|n`$ are excited states of $`\widehat{H}`$, and $`E_0`$, $`E_n`$ are their respective eigenenergies. Although $`\eta =0^+`$ is infinitesimal, a finite value may be used to broaden the resonances at $`\mathrm{}\omega =E_nE_0`$ and to reduce finite-size effects. $`\widehat{J}`$ is the current operator

$$\widehat{J}=\frac{iaet}{\mathrm{}}\underset{l;\sigma }{}\left(\widehat{c}_{l,\sigma }^+\widehat{c}_{l+1,\sigma }\widehat{c}_{l+1,\sigma }^+\widehat{c}_{l,\sigma }\right),$$
(4)
where $`e`$ is the charge of an electron. I set $`a=e=\mathrm{}=1`$ throughout, and $`t=1`$ is used in figures showing the optical conductivity, which means that $`\sigma _1(\omega )`$ is given in units of $`e^2a/\mathrm{}`$ and $`\omega `$ in units of $`t/\mathrm{}`$.

We note that the current operator is invariant under the spin-flip transformation but antisymmetric under charge conjugation and spatial reflection. Therefore, if the ground state $`|0`$ belongs to the symmetry subspace $`A_g^+(A_g,P_c,P_s)`$, only excited states $`|n`$ belonging to the symmetry subspace $`B_u^{}(B_u,P_c,P_s)`$ contribute to the optical conductivity. According to selection rules, the matrix element $`0|\widehat{J}|n`$ vanishes if $`|n`$ belongs to another symmetry subspace. In this paper, the excitation energy $`E_nE_0`$ of the lowest eigenstate with a non-zero matrix element $`0|\widehat{J}|n`$ is called the optical gap $`E_{\text{opt}}`$.

## III Numerical methods

DMRG steve ; dmrgbook  is known to be a very accurate method for one-dimensional quantum systems with short-range interactions such as the extended Hubbard Hamiltonian (1). In this work I use three different DMRG techniques to calculate ground states, excited states, and dynamic response functions. All three techniques are based on the finite system DMRG algorithm.

First, the usual ground-state DMRG method is used to calculate the ground state $`|0`$ for a fixed number $`N_\sigma `$ of electrons of each spin $`\sigma `$. This method provides the ground-state energy $`E_0(N_{},N_{})`$ and allows us to calculate ground-state expectation values $`0|\widehat{O}|0`$ for various operators $`\widehat{O}`$, such as static correlation functions. The (Mott) gap $`E_M`$ in the single-particle density of states of a Mott insulator can also be obtained using this approach. At half filling ($`N_\sigma =N/2`$) the Mott gap is simply given by

$$E_M=2[E_0(N_{}+1,N_{})E_0(N_{},N_{})]$$
(5)
because of the charge-conjugation symmetry.

Second, a symmetrized DMRG ramasesha  technique is used to calculate the lowest eigenstates $`|n`$ in the $`B_u^{}`$ symmetry sector. This method yields not only the eigenenergies $`E_n`$ of the lowest optically excited states (in particular, the optical gap $`E_{\text{opt}}`$), but also allow us to compute expectation values $`n|\widehat{O}|n`$ and thus to analyze the nature of these states. To optimize the DMRG program performance, my implementation of the charge-conjugation and spin-flip symmetries differs from the original idea presented in Ref. ramasesha, . This is explained in detail in appendix.

Finally, the dynamical DMRG method JGE ; jeckel  is used to compute the optical conductivity (2) convolved with a Lorentzian distribution of width $`\eta >0`$. Comparisons with exact results have shown that DDMRG is a very reliable numerical method, which yields spectra with an accuracy comparable to exact diagonalizations but for much larger systems. JGE ; EGJ ; jeckel

All DMRG methods have a truncation error which is reduced by increasing the number $`m`$ of density-matrix eigenstates kept (for more details, see Refs. steve,  and dmrgbook, ). Varying $`m`$ allows one to compute physical quantities (including spectra) for different truncation errors and thus to obtain error estimates on these quantities. I have systematically used this procedure to estimate the precision of my numerical calculations and adjusted the maximal number $`m`$ of density-matrix states to reach a desired accuracy. The largest number of density-matrix eigenstates used in this work is $`m=1000`$. For all numerical results presented in this paper DMRG truncation errors are negligible.

All numerical calculations have been performed on lattices with an even number $`N`$ of sites using open boundary conditions. As we are interested in the properties of the Hamiltonian (1) in the thermodynamic limit, numerical calculations have always been carried out for several system sizes $`N`$ in order to investigate finite-size effects. The largest system size used here is $`N=512`$. If necessary, the results have been extrapolated to the infinite system limit $`N\mathrm{}`$. To evaluate finite-size effects in a continuous spectrum one has to compute it for different sizes while keeping $`\eta N`$ = const. jeckel  In this work $`\eta N=12.8t`$ is used. For all numerical results presented in this paper finite-size effects (including chain-end effects) are negligible unless discussed explicitly. For spectra this means that finite-size effects are completely hidden by the broadening $`\eta `$. More precisely, DDMRG results for finite $`N=12.8t/\eta `$ are not distinguishable from the corresponding infinite-system spectra convolved with a Lorentzian distribution of width $`\eta `$ (see the discussion in Ref. jeckel, ).

## IV Results

To facilitate the comparison with analytical results, the discussion of optical excitations in the Mott insulating phase is divided in three subsections: the limit of a large Mott gap, the regime of finite Mott gaps, and the limit of a small Mott gap. Note, however, that the Mott gap $`E_M`$ just fixes the energy scale; the minimal energy required to create a charge excitation is $`E_M/2`$ but optical excitations do not differ qualitatively as $`E_M`$ varies if everything else is kept constant. In all cases, there is spin-charge separation and the spin sector is gapless. Elementary excitations in the charge sector are spinless bosons in the lower and upper Hubbard bands. Optical excitations are always made of an even number ($`2`$) of elementary excitations with opposite charges (to preserve charge neutrality). The different types of optical excitations and optical spectra found in the model (1) result from the residual interactions (essentially the non-local part of the Coulomb repulsion, here $`V`$) between the elementary charge excitations.

### IV.1 Limit of a large Mott gap

In the strong-coupling limit $`Ut`$, the properties of the model (1) in the Mott insulating phase can be described using simple concepts. In the ground state double occupation is prohibited and there is exactly one electron on each site. Elementary charge excitations can be represented as an empty site (holon in the lower Hubbard band) or a doubly occupied site (doublon in the upper Hubbard band). The minimal energy required to create a holon or a doublon is $`E_M/2=U/2O(t)t`$. Optical excitations always consist of an equal number of holons and doublons to conserve the total charge. The ionicity of excited states is defined as the change in the number of doubly occupied sites with respect to the ground state

$$I_n=n|\widehat{N}_d|n0|\widehat{N}_d|0,$$
(6)
where $`\widehat{N}_d=_l\widehat{n}_{l,}\widehat{n}_{l,}`$ and $`|0,|n`$ denote the ground state and excited states, respectively. Thus, $`I_n`$ is a measure of the number of doublons (or equivalently of holons) created by an excitation. Depending on the strength of the nearest-neighbor interaction parameter $`V`$, the low-energy optical excitations are made of a single doublon-holon pair ($`I_n=1`$) or are collective excitations of several such pairs ($`I_n>1`$). Note that $`I_n`$ is also equal to the derivative of the excitation energy $`E_nE_0`$ with respect to $`U`$.

#### Single holon-doublon pair

For $`V<U/3+O(t)`$, optical excitations consist of a single holon-doublon pair and the optical properties, which can be calculated exactly, florian ; gallagher  depend only on the parameters $`V`$ and $`t`$. For $`0V2t`$, holon and doublon are independent. A schematic representation of this state is shown in Fig. 1. This pair of free charge excitations gives rise to a continuous band in the optical spectrum $`\sigma _1(\omega )`$. The band starts at the Mott gap $`E_M=U4t`$ and has a width of $`8t`$. As there is no optical excitation with a lower energy than $`E_M`$, the Mott gap is also the optical gap $`E_{\text{opt}}`$. The optical spectra for $`V=0`$ and $`V=2t`$ are shown in Fig. 2 with a broadening $`\eta /t=0.1`$. At the conductivity threshold $`\sigma _1(\omega )`$ vanishes as $`\sqrt{\omega E_M}`$ for $`V<2t`$ but diverges as $`1/\sqrt{\omega E_M}`$ in the special case $`V=2t`$. The optical conductivity also has a small peak at $`\omega =UV`$ with 1 % of the spectral weight. EGJ  This peak is visible inside the band for $`V=0`$ and $`V=2t`$ in Fig. 2. It corresponds to a bound state made of dispersionless charge excitations florian  and can be seen as a localized exciton with size $`\xi =1`$.

For $`V>2t`$, there is also a continuous band starting at $`\omega =E_M=U4t`$ due to independent holon-doublon pairs, but the lowest optical excitation is now an exciton (i.e., a bound holon-doublon pair) with an energy $`\omega _{\text{exc}}=E_{\text{opt}}=UV4t^2/V`$ (Ref. florian, ). The term $`4t^2/V`$ is the kinetic energy lowering due to the exciton center-of-mass motion. Therefore, in the strong-coupling limit ($`Ut`$) the exciton binding energy is $`E_b=E_ME_{\text{opt}}=V4t+4t^2/V`$ ($`V>2t`$). This binding energy significantly differs from the incorrect result $`E_b=V`$ often reported in the literature, shuai  which is (approximately) valid only under the additional condition $`Vt`$. In the optical spectrum $`\sigma _1(\omega )`$ the exciton generates an isolated $`\delta `$-peak at $`\omega _{\text{exc}}`$ below the band onset. For $`Vt`$, an exciton is essentially the nearest-neighbor holon-doublon pair shown in Fig. 1. This is exactly the state generated by the current operator (4) applied to the ground state with one electron on each site. Thus the spectral weight is concentrated in the excitonic peak for $`Vt`$. This strong excitonic peak is already clearly visible for $`V=5t`$ in Fig. 2. For finite $`V/t`$, however, there is a finite probability of finding holon and doublon at a distance $`m>1`$. This probability can be calculated exactly EGJ

$$C(m)=C(1\delta _{m,0})e^{\kappa m}$$
(7)
with $`\kappa =2\mathrm{ln}(V/2t)`$ and a normalization constant $`C`$. The exciton size (the average holon-doublon distance) is then

$$\xi (V>2t)=\frac{V^2}{V^24t^2}$$
(8)
and decreases as $`V`$ increases. Correspondingly, one observes a progressive transfer of spectral weight from the band above $`E_M`$ to the excitonic peak at $`\omega _{\text{exc}}=UV4t^2/V`$ as $`V`$ increases. Note that, as for $`V2t`$, there is a small peak corresponding to a localized exciton at $`\omega =UV`$ in the optical conductivity. This peak lies in the band for $`V<4t`$ but is situated between the band and the strong excitonic peak for $`V>4t`$. In Fig. 2 it can be seen for $`V=5t`$ as a small bump at the foot of the strong peak.

Figure 2 shows optical spectra calculated with the DDMRG method on a 128-site lattice. On the scale of this figure, the DDMRG spectra are indistinguishable from the analytical results for an infinite system. florian ; EGJ  (For an expanded view showing small deviations, see Fig. 1 of Ref. EGJ, .) This perfect agreement confirms the validity of the strong-coupling calculations done in Ref. florian, . Moreover, it confirms once more that DDMRG can accurately reproduce infinite-system optical spectra. jeckel

#### Collective excitations

As we have just seen, the creation energy of an excited state with ionicity $`I_n=1`$ is either $`UO(t)`$ for an unbound holon-doublon pair or $`UVO(t^2/V)`$ for an exciton. Once a first excitation has been created, however, the creation of a second holon-doublon pair bound to the first excitation requires only an energy $`U2V`$. Therefore, when $`V`$ becomes large enough, the lowest optical excitations are bound states of $`n_{\text{exc}}`$ excitons, called excitonic strings. mazumdar ; pleutin  A biexciton ($`n_{\text{exc}}=2`$) is shown in Fig. 1 as an illustration (see also Ref. guo, ).

The excitation energy of a $`n_{\text{exc}}`$-exciton string is

$$E(n_{\text{exc}})=UV+(n_{\text{exc}}1)(U2V)O(t^2/V).$$
(9)
Here the correction of order $`t^2/V`$ corresponds to the kinetic energy lowering due to the center-of-mass motion. An $`n_{\text{exc}}`$-exciton string ($`n_{\text{exc}}2`$) appears in the low-energy excitation spectrum, around or below the onset $`E_M=U4t`$ of the band of free holons and doublons, if $`E(n_{\text{exc}})E_M`$ or

$$V\frac{n_{\text{exc}}1}{2n_{\text{exc}}1}U+O(t).$$
(10)
Thus, the biexciton becomes a low-energy excitation for $`VU/3`$ and longer excitonic strings ($`n_{\text{exc}}3`$) for larger $`V`$. The case $`n_{\text{exc}}=1`$ corresponds to the usual exciton, which appears in the low-energy spectrum as soon as $`V>2t`$ as discussed above. As an $`n_{\text{exc}}`$-exciton string is made of $`n_{\text{exc}}`$ doublons and holons bound together, it is a neutral excitation, its ionicity is $`I=n_{\text{exc}}`$, and its length is $`2n_{\text{exc}}1`$ in units of the lattice constant. Excitonic strings have been observed in the non-linear optical spectrum of quasi-one-dimensional neutral mixed-stack charge-transfer solids and are known to contribute to the non-linear optical conductivity of models such as the extended Hubbard model (1).guo ; mazumdar  Naively, one does not expect excitonic strings with $`I=n_{\text{exc}}2`$ to contribute to the linear optical spectrum $`\sigma _1(\omega )`$. In the limit $`U/tt`$ the current operator (4) creates at most one holon-doublon pair and thus in Eq. (3) the matrix elements $`n|\widehat{J}|0`$ between an excited state $`|n`$ and the ground state $`|0`$ must vanish if the ionicity (6) is larger than 1. Yet, we will see in the next section that excitonic strings with $`n_{\text{exc}}2`$ are visible in the linear optical conductivity of the extended Hubbard model for large but finite couplings $`U/t`$ and $`V/t`$. The reason is that for any finite $`t`$ there are quantum charge fluctuations (virtual holon-doublon pairs) in all eigenstates of the Hamiltonian (1) which leads to small but finite matrix elements $`n|\widehat{J}|0`$ even if the average ionicity $`I_n`$ of an excitation $`|n`$ exceeds one,

$$n|\widehat{J}|0\left(\frac{t}{U2V}\right)^{(I_n1)}.$$
(11)
Thus, for $`VU/3t`$ the low-energy optical spectrum consists of a strong excitonic peak at $`\omega =E_{\text{opt}}=UVO(t^2/V)`$ followed by several \[$`n_{\text{exc}}=2,3,\mathrm{}(UV)/(U2V)`$\] weaker isolated $`\delta `$-peaks with exponentially decreasing spectral weight at $`\omega E_{\text{opt}}+(n_{\text{exc}}1)(U2V)`$. All these peaks appears below (or about) the onset of a weak continuum due to free holons and doublons at $`\omega =E_M`$.

As long as $`U2Vt`$, excitonic strings retain a well-defined size represented by an integer number $`n_{\text{exc}}`$ because the kinetic energy lowering $`t`$ due to size fluctuations is much smaller than the energy cost $`(U2V)`$ per exciton in the string. Close to the phase boundary ($`U2V`$) between the CDW ground state and the Mott insulator, hirsch ; jeckel2  however, size fluctuations become advantageous. Thus, for $`U2Vt`$ low-energy excitations (of the Mott insulator) are CDW droplets, which can be understood as superpositions of excitonic strings of every size,

$$|\psi =c_1|n_{\text{exc}}=1+c_2|n_{\text{exc}}=2+c_3|n_{\text{exc}}=3+\mathrm{}$$
(12)
with a broad distribution of coefficients $`c_n`$ ($`_n|c_n|^2=1`$). For $`U2V0^+`$ the distribution becomes flat, i.e., $`|c_n|^2`$ const. (For comparison, an $`n_{\text{exc}}`$-exciton string can be described by the above state with $`|c_n|^21`$ for $`n=n_{\text{exc}}`$ and $`|c_n|^21`$ for $`nn_{\text{exc}}`$.) CDW droplets in the Mott insulating phase are the analogue to the SDW droplets in the CDW insulator discussed by Hirsch. hirsch  As excitonic strings, these CDW droplets are neutral excitations, but one can generalize the concept to CDW droplets carrying charges (see below). The average size $`r_{\text{CDW}}`$ of a CDW droplet is related to its ionicity by $`r_{\text{CDW}}=2I=2_ii|c_i|^2`$. Its excitation energy is

$$E(r_{\text{CDW}})=UV+\frac{r_{\text{CDW}}}{2}(U2V)st,$$
(13)
where $`s>0`$ and $`st`$ represents the kinetic energy lowering due to droplet size fluctuations and center-of-mass motion ($`s4`$ for $`U2V0^+`$). Contrary to excitonic strings, the ionicity $`I`$ of a CDW droplet is not a integer number but can take any value $`1`$. Therefore, for $`U2Vt`$ there is a band of CDW droplet excitations starting at $`UVst`$. Moreover, the matrix element $`n|\widehat{J}|0`$ for a CDW droplet $`|n`$ is essentially given by the overlap $`c_1`$ with the single exciton state in Eq. (12). Thus in this regime one expects that the CDW droplets give rise to a band in the optical spectrum $`\sigma _1(\omega )`$ starting at $`\omega =E_{\text{opt}}=UVst`$. This band lies below the Mott gap $`E_M`$. It should be noted that the Mott gap is determined by the excitation energy of unbound holons and doublons, $`E_MU4t`$, as long as $`U2Vt`$. For $`U2Vt`$, however, CDW droplets carrying a charge $`\pm e`$ \[a CDW droplet (12) bound to an extra holon or doublon\] have a lower energy than a bare holon or doublon and reduce the gap for charge excitations (5) to $`E_M=Us^{}t`$ with $`s^{}8`$ for $`VU/2`$.

In summary, in the limit of a large Mott gap ($`Ut`$) there are four distinct regimes corresponding to four types of excitations in the low-energy optical spectrum: independent charge excitations for $`V2t`$, excitons for $`V>2t`$ but $`VU/3`$, excitonic strings for $`VU/3t`$ but $`U2Vt`$, and CDW droplets for $`U2Vt`$. We will see that these excitations are also found in the optical spectrum of the extended Hubbard model for finite interaction strengths and Mott gaps.

### IV.2 Regime of finite Mott gaps

For finite coupling parameters $`U`$ and $`V`$ the low-energy optical properties of the extended Hubbard model (1) can be calculated using the ground-state and symmetrized DMRG methods presented in Sec. III. Figure 3 shows the Mott gap $`E_M`$ and the optical gap $`E_{\text{opt}}`$ as a function of $`V`$ for three values of $`U`$. Both gaps increase monotonically with $`U`$ but decrease with increasing nearest-neighbor interaction $`V`$shuai ; kancharla  For all values of $`U`$ the optical gap equals the Mott gap (in the thermodynamic limit) as long as $`V2t`$ but for larger $`V`$, $`E_{\text{opt}}`$ becomes smaller than $`E_M`$. This suggests that for all $`U>0`$ the low-energy excitations are unbound for $`V2t`$ and bound for $`V>2t`$ as in the $`Ut`$ limit. florian ; shuai  Obviously, the condition $`V>2t`$ can be realized only for relatively strong coupling ($`U4t`$) because the Mott insulating phase exists only for $`V`$ up to $`V_cU/2`$. The Mott gaps in Fig. 3 are initially almost constant as $`V`$ increases then diminish significantly close to the phase boundary $`V_c`$. This agrees with the strong-coupling analysis in the previous section, which suggests that $`E_M`$ is essentially independent of $`V`$ for $`V_cVt`$ but is reduced by quantity $`t`$ as $`V`$ approaches the critical value $`V_c`$. Note that on the critical line between the CDW and Mott insulating phases, both gaps vanish for $`U3t`$ while the Mott gap clearly remains finite for stronger coupling ($`U4t`$).

To determine the nature of the low-energy optical excitations I have calculated their ionicity (6). Figure 4 shows the ionicity $`I_1`$ of the first optically excited state (the $`1B_u^{}`$ state) as a function of $`V`$ for three values of $`U`$. In the half-filled Hubbard model ($`V=0`$) $`I_1`$ increases monotonically from 0 at $`U=0`$ to 1 for $`U\mathrm{}`$, reflecting the increasingly ionic nature of the elementary charge excitations in the lower and upper Hubbard bands. The ionicity increases slowly with $`V`$ and remains below or close to 1 for most couplings ($`U,V`$), which confirms that the corresponding optical excitations are made of a single pair of elementary charge excitations. In the regime $`U2V>4t`$, however, one observes a rapid but continuous increase of $`I_1`$ to values larger than 2 as $`VV_cU/2`$. This shows that the lowest optical excitation has become a CDW droplet. Looking at higher optical excitations, one finds the same qualitative behavior of the ionicity as a function of $`U`$ and $`V`$. Additionally, one observes the formation of excitonic string with integer $`I=n_{\text{exc}}2`$.

To determine whether the pair of elementary charge excitations form bound (exciton) or unbound states, one can calculate the average distance between both excitations using an exciton correlation function. joerg ; joerg2 ; EGJ  For $`V<2t`$ I have found that this average distance always diverges with increasing system size $`N`$. This result definitively confirms that in this regime an optical excitation is a pair of independent charge excitations, in agreement with the strong-coupling analysis. For $`V>2t`$, the average distance tends to a finite value for $`N\mathrm{}`$ as expected for an exciton. The exciton size $`\xi `$ determined with this procedure is shown in Fig. 5 as a function of $`V`$ for two finite values of $`U`$. The exciton size in the limit $`Ut`$, Eq. (8), is also plotted for comparison. The size $`\xi `$ increases and diverges as $`V`$ tends to $`2t`$, showing the unbinding of the exciton at $`V=2t`$. Note that for $`U=40t`$ the sizes measured with the exciton correlation function EGJ ; joerg ; joerg2  agree perfectly with Eq. (8).

One can gain some knowledge about the nature of optically excited states by looking at the scaling of their spectral weight

$$W_n=\frac{\pi }{Na}\frac{|n|\widehat{J}|0|^2}{E_nE_0}$$
(14)
\[see Eqs. (2) and (3)\] with the system size $`N`$. For $`V2t`$ I have found that the spectral weight of low-lying excitations vanishes for $`N\mathrm{}`$, which is expected for states belonging to a continuum. (As there is an infinite number of states in a continuous band, the spectral weight of each state must go to zero as $`N\mathrm{}`$, so that the total weight in any finite frequency interval remains finite.) For $`V>2t`$ (but outside the CDW droplet regime), I have found that the optical weight $`W_1`$ of the $`1B_u^{}`$ state tends to a finite value in the thermodynamic limit. This corresponds to a $`\delta `$-peak (with total weight $`W_1`$) in the optical spectrum as expected for an exciton or an excitonic string. In the CDW droplet regime, finite-size effects become large and complex and, in most cases, it has not been possible to determine the scaling of the spectral weights (14).

#### Optical spectra

The above analysis shows that low-energy optical excitations in the regime of finite Mott gaps are identical to those found in the limit of a large Mott gap and can be interpreted using the simple theory developed for the strong-coupling limit in Sec. IV.1. Turning next to the optical spectrum I have calculated $`\sigma _1(\omega )`$ for various parameters $`40tU3t`$ and $`U/2V0`$ using the DDMRG method. I have found that the optical spectra of systems with a finite Mott gap closely resemble those observed in the limit of a large Mott gap.

As a first example, the optical conductivity $`\sigma _1(\omega )`$ is shown in Fig. 6 for $`U=40t`$ and several values of $`V`$ representing the four different regimes: free charge excitations ($`V=0`$), excitons ($`V=5t`$), excitonic strings ($`V=16t`$), and CDW droplets ($`V=19.97t`$). For $`V=0`$ there is a single continuous band starting at $`\omega =E_{\text{opt}}=E_M=36.14t`$. At the band edges the optical conductivity vanishes as $`\sqrt{\omega E_M}`$ as discussed in Ref. JGE, .

For $`V=5t`$ a strong excitonic $`\delta `$-peak appears at $`\omega _{\text{exc}}=E_{\text{opt}}=34.39t`$ below the Mott gap $`E_M=36.13t`$. The exciton has a size $`\xi 1.2`$ in perfect agreement with Eq. (8). There is also a weak continuous band of free charge excitations above $`E_M`$, which is only visible as a high-frequency tail of the exciton peak in Fig. 6. The gap between the excitonic peak and the band is not visible in Fig. 6 because of the large broadening $`\eta /t=0.4`$ used here, but it can be checked with a scaling analysis for $`\eta 0`$ ($`N\mathrm{}`$) as discussed in Ref. jeckel, . The only qualitative difference between the present result for $`U=40t`$ and the corresponding result in the limit $`Ut`$ (see Fig. 1 in Ref. EGJ, ) is the absence of the weak peak associated with a localized exciton at $`\omega =UV`$. Nevertheless, this weak peak is not an artifact of the strong-coupling limit because its existence has been confirmed in the Hubbard model ($`V=0`$) down to $`U=4t`$ (Ref. JGE, ). The finite spectral weight carried by the localized exciton originates from a ground-state dimer-dimer correlation of the spin degrees of freedom. florian  In the strong-coupling limit ($`U2Vt`$) of the extended Hubbard model (1), the effective exchange coupling between nearest-neighbor spins depends on the occupation of the neighboring sites if $`V0`$. Thus the effective spin Hamiltonian is not the one-dimensional Heisenberg model, in general. Only for $`V=0`$ or in the limit $`UV`$, the effective spin Hamiltonian reduces to the Heisenberg model with a constant exchange coupling $`J=4t^2/U`$. In this case, the ground state has the relevant spin dimer-dimer correlations and the localized exciton carried a finite optical weight as explained in detail in Ref. florian, . For finite $`U`$ and $`V`$, however, the spin dimer-dimer correlation is presumably destroyed by the fluctuations of the spin exchange coupling and thus the optical weight of the localized exciton vanishes.

For $`V=16t`$ the condition $`VU/3`$ is satisfied and excitonic strings appear in the optical spectrum below or around the Mott gap $`E_M=35.58t`$. As seen in Fig. 6, most of the spectral weight is concentrated in the exciton of size $`\xi =1.0`$ at $`\omega _{\text{exc}}=E_{\text{opt}}=23.53tUV`$. The biexciton at $`\omega =32.35t2U3V`$ is barely visible in Fig. 6. The optical conductivity $`\sigma _1(\omega )`$ is again shown in Fig. 7(a) on a logarithmic scale. The isolated peaks associated with both excitations are now clearly visible. The measured ionicity (6) is $`I=1.1`$ and $`I=2.2`$ for the exciton and the biexciton, respectively. In Fig. 7(a) the remnant of the continuous band of free charge excitations and the triexciton (at $`\omega 3U5V=40t`$) are also visible in the interval $`\omega =3644t`$ above the Mott gap.

For $`V=19.97tU/2`$ the optical conductivity spectrum is radically different. The excitonic strings collapse into a band of CDW droplets with varying sizes. For instance, the $`1B_u^{}`$ state is a droplet of size $`r_{\text{CDW}}=8.9`$ with an energy $`E_{\text{opt}}15.5t`$. These CDW droplets give rise to a broad band in the optical conductivity spectrum shown in Fig. 6. The onset of this band is well below the Mott gap $`E_M=30.83t`$. On the logarithmic scale of Fig. 7(b), one sees that, in this particular case, the entire optical weight seems to be below $`E_M`$ (for $`\eta 0`$). The appearance of a band below the Mott gap is also visible in the current-current correlations for $`U=12t`$ and $`V=6t`$ presented in Ref. kancharla,  but the optical spectrum in the regime $`U2V`$ is not interpreted correctly in that work.

As a second example and to illustrate the finite-size-scaling analysis I have carried out for dynamical spectra, I discuss the optical conductivity $`\sigma _1(\omega )`$ for $`U=8t`$. Figure 8 shows the evolution of the optical conductivity for increasing nearest-neighbor repulsion $`V`$. For $`V=t`$ and $`V=2t`$, the spectrum contains a single continuous band due to free charge excitations starting at $`E_{\text{opt}}=E_M=4.67t`$ and $`4.53t`$, respectively. For $`V=3t`$ the spectrum consists of a strong peak corresponding to an exciton of size $`\xi =3.2`$ and energy $`\omega _{\text{exc}}=E_{\text{opt}}=3.86t`$, and of a weak band above the Mott gap $`E_M=4.10t`$. This band is visible in Fig. 8 as the high-frequency tail of the excitonic peak. For $`V=4t`$, CDW droplets of varying sizes dominate the optical spectrum. For instance, the $`1B_u^{}`$ state is a droplet of size $`r_{\text{CDW}}=5.6`$ with an excitation energy $`E_{\text{opt}}=1.55t`$ lower than the Mott gap $`E_M=2.29t`$. There is no intermediate regime with well-defined excitonic strings for this value of $`U`$.

The precise shape of $`\sigma _1(\omega )`$ cannot be determined from the sole results shown in Fig. 8 because of the finite resolution and system size used, $`\eta /t=12.8/N=0.1`$. To determine the properties of $`\sigma _1(\omega )`$ with maximal resolution ($`\eta 0`$) in the thermodynamic limit ($`N\mathrm{}`$), one can perform a scaling analysis with $`\eta N=const.`$ as explained in Ref. jeckel, . (Here I have used $`\eta N=12.8t`$.) The scaling analysis of the optical conductivity $`\sigma _1(\omega )`$ calculated with DDMRG always yields results which are qualitatively and quantitatively consistent with the properties of low-lying optical excitations determined using the ground-state and symmetrized DMRG methods. For instance, $`\sigma _1(\omega )`$ vanishes for all $`\omega <E_{\text{opt}}`$ and there is a continuous band for $`\omega E_{\text{opt}}`$ or a $`\delta `$-peak at $`\omega =E_{\text{opt}}`$ in the limit $`\eta 1/N0`$.

The scaling analysis of the conductivity maximum $`\sigma _{\text{max}}=\sigma _1(\omega _{\text{max}})`$ in DDMRG spectra is illustrated in Fig. 9 for the same interaction parameters as in Fig. 8. For $`V=t`$, $`\omega _{\text{max}}`$ tends to a value $`(5.15t)`$ larger than the optical gap $`E_{\text{opt}}=4.67t`$ calculated with symmetrized DMRG for $`N\mathrm{}`$ \[see Fig. 9(a)\] while $`\sigma _{\text{max}}`$ tends to a finite value for $`\eta 0`$ \[see Fig. 9(b)\]. Moreover, the derivative of $`\sigma _1(\omega )`$ has a maximum that diverges as $`1/\sqrt{\eta }`$ for $`\eta 0`$ at $`\omega =E_{\text{opt}}`$. These features correspond to a continuum that vanishes as $`\sqrt{\omega E_{\text{opt}}}`$ at the conductivity threshold and goes through a maximum just above the optical gap at $`\omega _{\text{max}}1.1E_{\text{opt}}`$ (see Ref. jeckel, ). For $`V=2t`$, $`\omega _{\text{max}}`$ tends to the same value as the optical gap $`E_{\text{opt}}`$ for $`N\mathrm{}`$ \[see Fig. 9(a)\] and $`\sigma _{\text{max}}`$ diverges as $`1/\sqrt{\eta }`$ for $`\eta 0`$ \[see Fig. 9(b)\]. These features correspond to a continuum that diverges as $`1/\sqrt{\omega E_{\text{opt}}}`$ at the conductivity threshold. jeckel  (Note that for all values of $`U`$ investigated $`\sigma _1(\omega )`$ displays this divergence at $`V=2t`$.) Therefore, the features of the optical spectrum for $`V2t`$ are similar to those found in the strong-coupling limit. For $`V=3t`$, $`\omega _{\text{max}}`$ tends for $`N\mathrm{}`$ to the same value as the optical gap $`E_{\text{opt}}`$ (which is smaller than the Mott gap in the thermodynamic limit) \[see Fig. 9(a)\] and $`\sigma _{\text{max}}`$ diverges as $`1/\eta `$ for $`\eta 0`$ \[see Fig. 9(b)\]. These features correspond to a $`\delta `$-peak at $`\omega =E_{\text{opt}}`$. Moreover, $`\sigma _1(\omega )`$ vanishes between $`E_{\text{opt}}`$ and $`E_M`$ but remains finite above $`E_M`$ in the limit $`\eta 1/N0`$. Therefore, in the thermodynamic limit the spectrum for $`V=3t`$ (shown in Fig. 8 for $`N=128`$ sites) consists of an excitonic $`\delta `$-peak separated from the band of independent charge excitations as in the strong-coupling limit. In the CDW droplet regime (i.e., close to the critical line $`V_cU/2`$ separating Mott and CDW phases) finite-size effects are more complicated and larger than in the other regimes. As a consequence, for $`V=4t`$ it has not been possible to perform a conclusive analysis with the largest system sizes ($`N=256`$) available. It seems that the low-energy spectrum contains a $`\delta `$-peak at $`E_{\text{opt}}`$ and a band starting immediately above $`E_{\text{opt}}`$, both due to CDW droplets.

In summary, I have found that the optical properties for all finite Mott gaps (i.e., for all $`U>0`$, $`V0`$ in the Mott phase) are qualitatively similar to those calculated in the limit of a large Mott gap (Sec. IV.1, Refs. EGJ, , florian, ):

(i) For $`V2t`$ (that is the only possible case for $`U4t`$), independent charge excitations give rise to a continuous band starting at the Mott gap $`E_M`$, which is equal to the optical gap. The band width is typically $`8t`$. For $`V<2t`$, $`\sigma _1(\omega )`$ vanishes smoothly at the threshold $`E_{\text{opt}}`$, typically as $`\sqrt{\omega E_{\text{opt}}}`$. At $`V=2t`$, $`\sigma _1(\omega )`$ diverges as $`1/\sqrt{\omega E_{\text{opt}}}`$ for $`\omega E_{\text{opt}}0^+`$.

(ii) For $`V>2t`$ but $`U>2V+O(t)`$ and $`V<U/3+O(t)`$ (this is possible only for $`U4t`$), the optical spectrum consists of an excitonic $`\delta `$-peak below the Mott gap and a band due to free charge excitations above $`E_M`$. Most of the optical weight is in the excitonic peak for $`V3t`$.

(iii) If $`U`$ is large enough ($`U12t`$) excitonic strings appear in the low-energy spectrum for $`V>U/3+O(t)`$ but $`U>2V+O(t)`$. They generate isolated $`\delta `$-peaks below the Mott gap $`E_M`$ in the optical conductivity $`\sigma _1(\omega )`$ with a separation between peaks of $`\mathrm{\Delta }\omega U2V`$. The first peak is an exciton and contains most of the spectral weight. A very weak band due to free charge excitations still exists above $`E_M`$.

(iv) Close to the boundary $`V_cU/2`$ between the Mott and CDW phases, if $`V`$ exceeds $`2t`$, CDW droplets dominate the low-energy spectrum and give rise to a broad band (including sharp peaks) starting below the gap for charge excitations (5).

### IV.3 Limit of a small Mott gap

In the limit of a small Mott gap ($`E_Mt`$) the coherence length $`4t/E_M`$ becomes very large and it is not possible to carry out numerical simulations on lattices large enough ($`N4t/E_M`$) to determine the optical spectrum with confidence. Fortunately, in this limit field-theoretical methods provide generic results for the low-energy optical spectrum of a one-dimensional Mott insulator. JGE ; controzzi ; EGJ  Field-theoretical results are applicable to lattice models such as Eq. (1) for gaps up to $`E_Mt`$, which makes possible a direct quantitative comparison of field theory and DDMRG calculations. JGE ; EGJ ; jeckel

In the field-theoretical approach, elementary charge excitations are holons (in the lower Hubbard band) and anti-holons (in the upper Hubbard band). Optical excitations are made of a equal number of holons and anti-holons. Assuming that the low-energy excitations consist of one holon-antiholon pair, the optical conductivity is

$$\sigma _1^{\text{FT}}(\omega )=AS_\beta (\omega /E_M)$$
(15)
where $`S_\beta (x)`$ is a known function depending on the field-theory interaction parameter $`0<\beta 1`$, and $`A`$ is a unknown constant which sets the conductivity scale. Strictly speaking, this result is exact only for $`\omega <2E_M`$ and $`\beta ^2>1/3`$ but it has been found by comparison with DMRG results that corrections for $`\omega >2E_M`$ are usually negligible. JGE ; EGJ  For $`\beta ^21/2`$, $`\sigma _1^{\text{FT}}(\omega )`$ describes a single continuous band starting at $`E_{\text{opt}}=E_M`$, which is due to independent holons and anti-holons. The optical conductivity vanishes smoothly as $`\sqrt{\omega E_{\text{opt}}}`$ for $`\beta ^2>1/2`$ and diverges as $`1/\sqrt{\omega E_{\text{opt}}}`$ for $`\beta ^2=1/2`$ at the conductivity threshold. For $`1/3<\beta ^2<1/2`$, there is a $`\delta `$-peak at $`\omega =E_{\text{opt}}<E_M`$ in addition of the band starting at $`\omega =E_M`$. The $`\delta `$-peak is due to a bound holon-antiholon pair (exciton). For $`\beta ^2<1/3`$ additional excitons and excitonic strings (made of several holon-antiholon pairs) appear in the spectrum. Therefore, field-theoretical predictions for the optical conductivity of a one-dimensional Mott insulator are qualitatively similar to what we have found in the extended Hubbard model (1) using a strong-coupling analysis and DDMRG simulations.

The field-theory parameters $`E_M`$, $`\beta `$, and $`A`$ must be estimated numerically by comparison with DMRG results because one does not know their relations to the lattice model parameters $`U,V,t`$. I have first determined several couplings $`(U,V)`$ which yield approximately the same Mott gap $`E_M/t0.60.7`$. These couplings are listed in Table 1 with the corresponding Mott and optical gaps calculated using the ground-state and symmetrized DMRG methods. Then I have calculated the optical conductivity $`\sigma _1(\omega )`$ for these parameters using DDMRG. Some results are shown in Fig. 10. Note the progressive displacement of spectral weight to lower energy as $`V`$ increases although the Mott gap remains almost constant (see Table 1). To determine the parameters $`\beta `$ and $`A`$ one can now compare the field-theoretical spectra with the DDMRG data. \[$`\sigma _1^{\text{FT}}(\omega )`$ has to be convolved with a Lorentzian distribution of appropriate width $`\eta `$ to make a direct comparison. jeckel \] This procedure yields the parameter $`\beta ^2`$ listed in Table 1. As expected the boundaries between independent charge excitations and excitons in the field theory and in the lattice model (1) coincide: $`\beta ^2=1/2`$ corresponds to $`V=2t`$.

For $`V2t`$ (corresponding to $`\beta ^21/2`$) one can find parameters so that $`\sigma _1^{\text{FT}}(\omega )`$ perfectly fits the numerical data over a wide frequency range. For instance, in Fig. 11 no difference is visible between the DDMRG spectrum for $`U=4.15t`$ and $`V=2t`$ and the fitted field-theoretical spectrum up to $`\omega =2t3E_M`$. For $`V>2t`$, however, discrepancies between DDMRG and field-theory results appear and grow progressively stronger as $`V`$ increases. It is no longer possible to find parameters $`\beta `$ and $`A`$ to reproduce the DDMRG spectra over a significant frequency range above $`\omega =E_{\text{opt}}`$. Instead $`\beta `$ and $`A`$ are set by the optical gap and the total spectral weight. This yields the values of $`\beta `$ shown in parenthesis in Table 1. As an example, one see in Fig. 11 that the field-theoretical spectrum differs significantly from the DDMRG result for $`U=6t`$ and $`V=3.115t`$, although optical gap, Mott gap, and total spectral weight are identical for both spectra. The DDMRG result shows that there is substantial optical weight both at $`E_{\text{opt}}=0.445t`$ and above the Mott gap $`E_M=0.643t`$ while, according to field theory, EGJ  for a ratio $`E_{\text{opt}}/E_M0.7`$ the optical conductivity should be dominated by an excitonic peak at $`\omega =E_{\text{opt}}`$ with very little weight in the holon-antiholon band above the Mott gap $`E_M`$ and in the exciton-exciton continuum above $`\omega 1.7E_{\text{opt}}1.1t`$.

This disagreement between field theory and DDMRG results is not due to a failure of either method in the excitonic regime. It has been shown that field theory and DDMRG calculations for excitons agree very well in the extended Hubbard model with next-nearest-neighbor repulsion. EGJ  The problem is that the field-theory approach assumes that the low-energy optical excitations are made of two elementary charge excitations. In the extended Hubbard model (1), however, the conditions $`E_M4t`$ and $`V>2t`$ are satisfied only close to the phase boundary $`V_c=U/2`$, where low-energy excitations are CDW droplets made of many elementary charge excitations. For instance, the ionicity of the $`1B_u^{}`$ state becomes significantly larger than 1 as $`V`$ increases above $`2t`$ as seen in Table 1. Therefore, the field-theory approach is not applicable to the lattice model (1) with $`V>2t`$ even in the limit of a small Mott gap. Moreover, the extended Hubbard model (1) cannot describe a Mott insulator with a small gap and an exciton in the optical spectrum for any parameters $`U`$ and $`V`$.

## V Conclusion

I have investigated the linear (one-photon) optical excitations of a one-dimensional Mott insulator, the half-filled extended Hubbard model, using DMRG methods. Four types of optically excited states have been found: pairs of free (unbound) charge excitation, excitons, excitonic strings, and CDW droplets. Correspondingly, there are four different regimes in the model parameter space $`(U,V)`$ depending on the nature of the low-energy optical excitations. They are shown in the schematic “phase diagram” of Fig. 12. Note that only the $`V=2t`$ line separating the regime of free excitations from that of bound excitations represents a sharp transition. The other dashed lines represent smooth crossover from one regime to another. In each regime one observes optical spectra with distinct features. In all cases, optical excitations are made of an even number of elementary excitations carrying opposite charges in the lower and upper Hubbard bands. The different types of excitations and optical spectra found in this system result from the residual interactions between these elementary excitations.

For $`V2t`$ the low-energy optical excitations are made of two unbound elementary charge excitations. They give rise to a single continuous band in the optical spectrum, starting at the Mott gap $`E_M`$, which thus equals the optical gap $`E_{\text{opt}}`$. The optical conductivity vanishes smoothly as $`\sqrt{\omega E_{\text{opt}}}`$ at the threshold $`E_{\text{opt}}`$, except for $`V=2t`$, where it diverges as $`1/\sqrt{\omega E_{\text{opt}}}`$. For $`V>2t(U4t)`$ but $`V<U/3+O(t)`$ and $`U2Vt`$, the lowest optical excitation is an exciton (a neutral excitation made of two bound elementary charge excitations) at an energy $`E_{\text{opt}}`$ lower than the Mott gap $`E_M`$. This exciton gives rise to an isolated $`\delta `$-peak in the optical conductivity at $`\omega =E_{\text{opt}}`$. Furthermore, one still finds a continuous band starting at the Mott gap $`E_M`$ due to free charge excitations. For $`V>U/3+O(t)`$ and $`V>2t`$ but $`U2Vt`$ ($`U12t`$) the low-energy optical excitations are excitonic strings consisting of $`n_{\text{exc}}1`$ excitons with energy $`E(n_{\text{exc}})E_{\text{opt}}+(n_{\text{exc}}1)(U2V)`$. The spectrum consists of several isolated $`\delta `$-peaks at $`\omega =E(n_{\text{exc}})<E_M`$. In this regime one still observes a very weak band of free charge excitations starting at $`E_M`$ in the spectrum. Finally, close to the CDW phase boundary ($`U2Vt`$) for $`V>2t`$ the low-energy excitations are CDW droplets and give rise to a broad band starting below the Mott gap $`E_M`$.

As long as optical excitations are made of a pair of (bound or unbound) elementary charge excitations (i.e, the excitation ionicity is $`I1`$), the optical spectra calculated numerically with DDMRG agree perfectly with the analytical results obtained with a strong-coupling analysis or with field-theoretical methods. This agreement confirms the accuracy and the power of the DDMRG method for calculating dynamical spectra in the thermodynamic limit. It also confirms the wide range of validity of both analytical approaches.

Some results presented here suggest further investigations. First, in the extended Hubbard model (1) no exciton exists in the regime of a small Mott gap ($`E_Mt`$), which is relevant for some real materials such as conjugated polymers. kiess ; sari  It is believed that an electron-electron interaction with a longer range EGJ ; gallagher ; barford  or a lattice dimerization shuai ; joerg2  can lead to the formation of excitons in systems with small gaps. However, the precise nature of the optical excitations in such systems is still controversial. robert  The approach used here for the extended Hubbard model will enable us to determine the optical properties of these systems reliably. Second, excitonic strings appear in the linear optical conductivity spectrum for strong interaction $`U,Vt`$ because of the weak hybridization of excitonic strings with different sizes $`n_{\text{exc}}`$. Experimentally, excitonic strings have been observed in the non-linear optical absorption only. mazumdar  It would be desirable to check if excitonic strings can be found in the linear optical absorption of materials which are believed to be large-gap one-dimensional Mott insulators such as Cu oxides and Ni halides. hasan ; kishida ; mizuno  Last, there is a clear boundary between free and bound excitations in the low-energy optical spectrum at $`V=2t`$. As discussed in Ref. jeckel2,  the nature of the low-energy charge excitation seems to be correlated with the order of the transition from the Mott insulating phase to the CDW insulating phase. It is likely that the tricritical point where the transition changes from continuous to first order is located precisely on the line $`V=2t`$. This suggests the existence of a hidden symmetry in the charge sector of the extended Hubbard model (1) at $`V=2t`$. It would be interesting to investigate this feature further.

###### Acknowledgements.

I gratefully acknowledge helpful discussions with R. Bursill, F. Essler, and S. Pleutin and I thank F. Gebhard for his support and many stimulating conversations.

*

## Appendix A

In the original implementation of the charge-conjugation and spin-flip symmetries for DMRG calculations, ramasesha  an explicit matrix representation of the superblock Hamiltonian is built. This matrix can be projected onto a symmetry subspace with chosen parities $`P_c`$ and $`P_s`$, which allows one to compute eigenstates of this symmetry and reduces the computer memory and CPU time required. In an efficient implementation of DMRG, however, an explicit representation of the superblock Hamiltonian should not be constructed (at least for quasi-one-dimensional systems with only short-range interactions). A representation in terms of tensor products of matrices uses much less memory and is also much faster. reinhard  Projecting this representation onto a symmetry subspace slows down the program considerably. Therefore, instead of a projection, I use an exact diagonalization technique xenophon  to shift the chosen symmetry subspace to lower energy.

Let $`\widehat{P}_c`$ and $`\widehat{P}_s`$ be the charge-conjugation and spin-flip operators for the full lattice with eigenvalues $`P_c=\pm 1`$ and $`P_s=\pm 1`$. As $`\widehat{P}_c`$ and $`\widehat{P}_s`$ commute with the Hamiltonian $`\widehat{H}`$, the operator

$$\widehat{H}^{}=\widehat{H}\lambda _c\widehat{P}_c\lambda _s\widehat{P}_s$$
(16)
has the same eigenstates as $`\widehat{H}`$ but its eigenvalues are shifted, $`E_n^{}=E_n\pm \lambda _c\pm \lambda _s`$, where the signs $`\pm `$ are given by the eigenstate parities $`P_c`$ and $`P_s`$. It is obvious that the lowest eigenstates of $`\widehat{H}^{}`$ lies in the symmetry subspace with $`P_c=\lambda _c/|\lambda _c|`$ and $`P_s=\lambda _s/|\lambda _s|`$ provided $`|\lambda _c|`$ and $`|\lambda _s|`$ are large enough. Therefore, one can simply apply the usual ground-state DMRG approach to the Hamiltonian $`\widehat{H}^{}`$ with appropriate values of $`\lambda _c`$ and $`\lambda _s`$ to obtain the lowest eigenstates in any symmetry sector. A similar approach has already been used to shift states with high total spin $`S`$ to higher energy in a DMRG calculation. stephane  Using the method proposed recently for including a non-abelian symmetry group in a DMRG calculation would be a further improvement. ian