text:
# Force correlations in the 𝑞–model for general 𝑞–distributions

## I Introduction

One of the main challenges of granular media is to characterize the network of microscopic forces in a static bead pack. In order to describe the corresponding force fluctuations, Liu et al. liu  introduced the $`q`$–model. In this model, the beads are placed on a regular lattice and the (scalar) forces are stochastically transmitted, by random fractions denoted by the symbol $`q`$. Even in its simplest version, where one assumes a uniform $`q`$–distribution, it already reproduces the main feature of the experimental observations: the probability for large forces decays exponentially liu ; exp1 ; exp2 . Although for this uniform $`q`$–distribution the forces become totally uncorrelated, in general, correlations do persist cop . In the present study, we investigate for which $`q`$–distributions this is the case and we reveal the surprising nature of these correlations. In order to perform an analytical study, we restrict ourselves to the scalar $`q`$–model and allow only correlations between $`q`$–values in a layer. More sophisticated lattice models, which include the vector nature of the force and allow correlations between layers are not considered here vec .

Although the $`q`$–model is particularly simple, its behavior turns out to be very rich. First of all, there is a so-called critical $`q`$–distribution, which produces a force distribution that decays algebraically instead of exponentially cop ; bouch . It therefore forms a critical point in the space of $`q`$–distributions, and its properties were recently investigated in great detail raj ; lew . A second intriguing issue concerns the top–down dynamics of force correlations (the downward direction can be interpreted as time) raj ; lew ; jac . Even if both in the initial state (top layer) and in the asymptotic state (infinite depth) all forces are uncorrelated, there will be correlations at all intermediate levels. Correlations become longer in range while their amplitudes diminish in a diffusion process, and as a result, the asymptotic force distribution is only approached algebraically jac . This process is closely related to the subject of this study, namely the presence of force correlations at infinite depth.

Let us recapitulate the definition of the $`q`$–model. The beads are assumed to be positioned on a regular lattice. Let $`f_i`$ be the force in the downward direction on the $`i`$th bead in a layer. This bead makes contact with a number of $`z`$ beads in the layer below, which we indicate by the indices $`i+\alpha `$. The $`\alpha `$’s are displacement vectors in the lower layer as shown in Fig. 1. Bead $`i`$ transmits a fraction $`q_{i,\alpha }`$ of the force $`f_i`$ to the bead $`i+\alpha `$ underneath it. These fractions are taken stochastically from a distribution satisfying the constraint

$$\underset{\alpha }{}q_{i,\alpha }=1,$$
(1)
which assures mechanical equilibrium in the vertical direction. So, we can write the force $`f_j^{}`$ on the $`j`$th bead in a layer as

$$f_j^{}=\underset{\alpha }{}q_{j\alpha ,\alpha }f_{j\alpha }.$$
(2)
As the weights of the particles are unimportant at infinite depth, we have left out the so–called injection term. The distribution of forces at infinite depth depends on the $`q`$–distribution $`H(\stackrel{}{q})`$, where the symbol $`\stackrel{}{q}`$ is a shorthand for all the $`q_{i,\alpha }`$ at a given layer. This $`H(\stackrel{}{q})`$ can be any function that is constrained by Eq.(1). If we now assume that there are no correlations between the $`q`$–values at different sites, the $`q`$–distribution is of the form

$$H(\stackrel{}{q})=\underset{i}{}\eta (\stackrel{}{q}_i)\delta \left(1\underset{\alpha }{}q_{i,\alpha }\right),\stackrel{}{q}_i=\{q_{i,\alpha }\},$$
(3)
where $`\eta (\stackrel{}{q}_i)`$ is symmetric in its arguments $`q_{i,\alpha }`$. Although we will refer to these $`q`$–distributions as “uncorrelated”, note that there are always correlations between the $`q_{i,\alpha }`$ of the same site due to the $`\delta `$–constraint.

In the first part of this study, we show that there is only a limited set of $`\eta (\stackrel{}{q}_i)`$ for which the stationary force distribution can be written as a product of single–site distributions, and therefore is totally uncorrelated. This set is an extension of the set that was already identified by Coppersmith et al. cop . In their extensive study, they also provided numerical evidence that, in general, correlations can persist. We will show that correlations are still absent in the second order moments. However, higher order correlations do exist and surprisingly enough, these turn out to decay algebraically. The results for the triangular packing and the fcc packing are summarized in Table 1, section VII. In the last part of this work, we show that one induces 2–point force correlations by allowing correlations between $`q`$–values on different sites in a layer. These correlations will generically vanish with a power law, except for the triangular packing, where the decay of force correlations follows the decay of the $`q`$–correlations.

The paper is organized as follows. In section II we derive a criterion that a distribution $`\eta (\stackrel{}{q}_i)`$ has to obey in order to produce an uncorrelated stationary state. We then show in section III, that this criterion is only obeyed for a limited set of $`\eta (\stackrel{}{q_i})`$. After that, we study the nature of the correlations, by writing the evolution of the force moments as master equations in section IV, and by analysing the stationary solutions of these equations in section V. Section VI deals with the effects of allowing correlations between the $`\stackrel{}{q_i}`$ of different sites in a layer, and the paper closes with a discussion.

## II Criterion for factorization

Using the recursive nature of the force transmission, Eq.(2), one can write down the following recursive relation for the force distribution cop ; jac :

$`P^{}(\stackrel{}{f}^{})`$ $`=`$ $`{\displaystyle H(\stackrel{}{q})𝑑\stackrel{}{q}P(\stackrel{}{f})𝑑\stackrel{}{f}}`$ (4)
$`\times {\displaystyle \underset{j}{}}\delta (f_j^{}{\displaystyle \underset{\alpha }{}}q_{j\alpha ,\alpha }f_{j\alpha }),`$
where we have introduced a vector notation for the forces in one layer $`\stackrel{}{f}=(f_1,\mathrm{},f_N)`$, and for the integrations we use the abbreviations

$$𝑑\stackrel{}{f}=\underset{i}{}_0^{\mathrm{}}𝑑f_i,$$
(5)
$$𝑑\stackrel{}{q}=\underset{i}{}𝑑\stackrel{}{q}_i=\underset{i}{}\underset{\alpha }{}_0^1𝑑q_{i,\alpha }.$$
(6)
It is often convenient to work with the Laplace transform of Eq.(4). Defining the Laplace transform as

$$\stackrel{~}{P}(\stackrel{}{s})=𝑑\stackrel{}{f}\mathrm{exp}(\stackrel{}{s}\stackrel{}{f})P(\stackrel{}{f}),$$
(7)
the recursion simplifies to cop ; jac

$$\stackrel{~}{P}^{}(\stackrel{}{s})=H(\stackrel{}{q})𝑑\stackrel{}{q}\stackrel{~}{P}(\stackrel{}{s}(\stackrel{}{q})),$$
(8)
with

$$s_i(\stackrel{}{q})=\underset{\alpha }{}q_{i,\alpha }s_{i+\alpha }.$$
(9)
The two representations Eq.(4) and (8) are equivalent, and they will both be used, depending on the nature of the problem.

The force distribution at infinite depth $`P^{}(\stackrel{}{f})`$ or $`\stackrel{~}{P}^{}(\stackrel{}{s})`$ can be obtained by finding the fixed point of the recursive relation. The main question of this section is to determine whether a given $`H(\stackrel{}{q})`$ leads to a $`P^{}(\stackrel{}{f})`$ that is simply a product of single–site force distributions $`p^{}(f_i)`$. In section VI we will show that this can only be the case for $`q`$–distributions of the type Eq.(3). So for this section, the question is: which $`\eta (\stackrel{}{q_i})`$ lead to uncorrelated asymptotic states?

To answer this question, let us assume that such a fixed point exists, i.e.

$$P^{}(\stackrel{}{f})=\underset{i}{}p^{}(f_i),\mathrm{or}\stackrel{~}{P}^{}(\stackrel{}{s})=\underset{i}{}\stackrel{~}{p}^{}(s_i).$$
(10)
Inserting this Ansatz into the Laplace representation of the recursion relation, Eq.(8), yields

$`\stackrel{~}{P}^{}(\stackrel{}{s})`$ $`=`$
$`{\displaystyle \underset{i}{}}{\displaystyle }`$ $`\eta (\stackrel{}{q}_i)`$ $`\delta \left(1{\displaystyle \underset{\alpha }{}}q_{i,\alpha }\right)d\stackrel{}{q}_i\stackrel{~}{p}^{}\left({\displaystyle \underset{\alpha }{}}q_{i,\alpha }s_{i+\alpha }\right)`$ (11)
$`=`$ $`{\displaystyle \underset{i}{}}\stackrel{~}{\psi }(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z}),`$
where the function $`\stackrel{~}{\psi }(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z})`$ is the outcome of the integral over the $`\stackrel{}{q_i}`$. The arguments represent the $`z`$ sites that are connected to site $`i`$ in the previous layer. Integrating out all forces except those at the $`z`$ sites connected to $`i`$ means putting all $`s_j=0`$ except the set $`\{s_{i+\alpha }\}`$:

$`\stackrel{~}{P}^{}(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z})`$ $`=`$ $`\stackrel{~}{\psi }(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z})`$ (12)
$`\times {\displaystyle \underset{\alpha }{}}\stackrel{~}{\psi }(s_{i+\alpha },0,\mathrm{})^{z1}.`$
This projection of the total force distribution can only factorize if $`\stackrel{~}{\psi }(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z})`$ is a product function as well, i.e.

$$\stackrel{~}{\psi }(s_{i+\alpha _1},\mathrm{},s_{i+\alpha _z})=\underset{\alpha }{}\stackrel{~}{\psi }(s_{i+\alpha }).$$
(13)
This leads to the following criterion for asymptotic factorization:

* Given a $`q`$–distribution $`\eta (\stackrel{}{q})`$, one can construct a factorized fixed point if, and only if, there is a function $`\stackrel{~}{\psi }(s)`$ that satisfies the following condition:

$$\eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}\left[\stackrel{~}{\psi }\left(\underset{\alpha }{}q_\alpha s_\alpha \right)\right]^z=\underset{\alpha }{}\stackrel{~}{\psi }(s_\alpha ).$$
(14)
* This function $`\stackrel{~}{\psi }(s)`$ is related to the single–site distribution as

$$\stackrel{~}{p}^{}(s)=\left[\stackrel{~}{\psi }(s)\right]^z.$$
(15)

Here, we omitted the site index $`i`$, and furthermore, our formulation depends only on $`z`$ (the number of $`q`$–values per site) and not on the details of the lattice.

## III Special class of $`q`$–distributions leading to factorization

It is a well–known fact that the so–called uniform distribution, in which $`\eta (\stackrel{}{q}_i)`$ is a constant, produces an uncorrelated asymptotic force distribution. In fact, Coppersmith et al. identified a countable set of $`q`$–distributions, of which the uniform distribution is a member, that have this property cop . Although it might seem obvious that a uniform distribution leads to an uncorrelated asymptotic state, it is really not trivial. Due to the constraint of Eq.(1), there are correlations between the $`q_{i,\alpha }`$ on each site $`i`$, which induce force correlations that only disappear under the special conditions discussed in the previous section, Eq.(14). In this section, we will show when these special conditions are obeyed.

There is a mathematical relation that is extremely important for the $`q`$–model zinn :

$`{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{1}{(1+s_\alpha )^r}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(zr)}{[\mathrm{\Gamma }(r)]^z}}{\displaystyle 𝑑\stackrel{}{q}\delta \left(1\underset{\alpha }{}q_\alpha \right)}`$
$`\times {\displaystyle \underset{\alpha }{}}(q_\alpha )^{r1}{\displaystyle \frac{1}{(1+_\alpha q_\alpha s_\alpha )^{zr}}}.`$
It holds for any real $`r>0`$. From this relation, it is immediately clear that for all $`q`$–distributions of the type

$$\eta (\stackrel{}{q})=\frac{\mathrm{\Gamma }(zr)}{[\mathrm{\Gamma }(r)]^z}\underset{\alpha }{}(q_\alpha )^{r1},r>0$$
(17)
there is a $`\stackrel{~}{\psi }(s)`$ that obeys Eq.(14), namely

$$\stackrel{~}{\psi }(s)=\frac{1}{(1+s)^r}.$$
(18)
The corresponding single–site force distributions are

$$\stackrel{~}{p}^{}(s)=\frac{1}{(1+s/zr)^{zr}}\mathrm{or}p^{}(f)=\frac{(zr)^{zr}}{\mathrm{\Gamma }(zr)}f^{zr1}e^{zrf}.$$
(19)
We rescaled the Laplace variable $`s`$, in order to put $`f=1`$. Coppersmith et al. already found these $`q`$–distributions for integer values of $`r`$, also based on Eq.(III) cop . However, it holds for any real $`r>0`$. This means that the set for which the stationary force distribution factorizes is substantially larger; it ranges from the infinitely sharp distribution ($`r\mathrm{}`$) to the critical distribution ($`r0`$) foot3 . Note that one recovers the results for the uniform distribution by putting $`r=1`$.

Although there is a huge variety of $`q`$–distributions that lead to uncorrelated force distributions, in general one cannot find a $`\stackrel{~}{\psi }(s)`$ that obeys Eq.(14). We will prove this by making a Taylor expansion of $`\stackrel{~}{\psi }(s)`$

$$\stackrel{~}{\psi }(s)=\underset{n=0}{\overset{\mathrm{}}{}}\psi _ns^n,$$
(20)
and then try to solve for the coefficients $`\psi _n`$ by imposing Eq.(14). It turns out that the equations can only be solved under special conditions, which are precisely obeyed by the class of $`q`$–distributions given by Eq.(17).

Let us first focus on the left hand side (LHS) of Eq.(14). The Taylor expansion will give rise to terms of the type $`(q_1s_1)^{n_1}(q_2s_2)^{n_2}\mathrm{}(q_zs_z)^{n_z}`$, which have to be integrated over all $`q_\alpha `$. This leads to terms $`s_1^{n_1}s_2^{n_2}\mathrm{}s_z^{n_z}`$ with prefactors given by the moments of $`\eta (\stackrel{}{q})`$

$$\overline{q_1^{n_1}q_2^{n_2}\mathrm{}q_z^{n_z}}=\eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}q_1^{n_1}q_2^{n_2}\mathrm{}q_z^{n_z}.$$
(21)
These moments are not independent, due to the constraint Eq.(1). In appendix A, we show that the moments

$$\eta _n=\eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}q_1^n$$
(22)
are in fact sufficient to characterize all relevant moments of Eq.(21). Besides the moments, there are of course additional prefactors consisting of combinations of the $`\psi _n`$; these are the quantities we try to find, for a given $`q`$–distribution $`\eta (\stackrel{}{q})`$.

The right hand side (RHS) of Eq.(14) also produces terms $`s_1^{n_1}s_2^{n_2}\mathrm{}s_z^{n_z}`$, with prefactors $`\psi _{n_1}\psi _{n_2}\mathrm{}\psi _{n_z}`$. The remaining task is to equate the prefactors of the terms $`s_1^{n_1}s_2^{n_2}\mathrm{}s_z^{n_z}`$ on both sides of the equation. This gives a set of equations, from which one can try to solve for the $`\psi _n`$.

The zeroth order equation is trivially obeyed for any $`\psi _0`$, as can be seen by putting all $`s_\alpha =0`$. For convenience we fix $`\psi _0=1`$. The same happens at first order, since for each $`\alpha `$, the LHS contains $`z`$ terms $`\psi _1\overline{q_\alpha }s_\alpha =1/z\psi _1s_\alpha `$, and the RHS is simply $`\psi _1s_\alpha `$. The first non-trivial equation appears at second order. There are two equations, for $`s_\alpha ^2`$ and for $`s_\alpha s_\alpha ^{}`$ where $`\alpha \alpha ^{}`$:

$$\{\begin{array}{ccc}\hfill \left(z\psi _2+\frac{z(z1)}{2}\psi _1^2\right)\eta _2& =& \psi _2\hfill \\ \hfill \left(z\psi _2+\frac{z(z1)}{2}\psi _1^2\right)\frac{2(1z\eta _2)}{z(z1)}& =& \psi _1^2\hfill \end{array}$$
(23)
Due to the constraint $`_\alpha q_{i,\alpha }=1`$, one can obtain an identity by multiplying the first equation by $`z`$, and adding it to the second equation multiplied by $`z(z1)/2`$. Hence, the two equations are not independent and $`\psi _2`$ can be solved. The value of $`\psi _2`$ depends only on $`\eta _2`$, the second moment of the $`q`$–distribution foot1 .

Working out the combinatorics of the higher orders, one finds the following general mathematical structure:

* At the $`n`$th order, there are as many equations as there are different partitions $`\{n_1,n_2,\mathrm{},n_z\}`$ that make $`_\alpha n_\alpha =n`$. Permutations should not be considered as different because $`\eta (\stackrel{}{q})`$ is symmetric in its arguments.
* One of these equations is dependent, as one can obtain an identity by adding the equations, after multiplication by appropriate factors.

For $`z=2`$, there are two third order equations, corresponding to the partitions $`\{3,0\}`$ and $`\{2,1\}`$, of which only one is independent. This means that $`\psi _3`$ can be solved as a function of $`\eta _2`$ (in appendix A we show that $`\eta _3`$ depends on $`\eta _2`$, for $`z=2`$). We run into problems at fourth order, where we have $`\{4,0\}`$, $`\{3,1\}`$ and $`\{2,2\}`$, and hence two a priori independent equations for one coefficient $`\psi _4`$. It turns out that the remaining equations are only identical if there is a relation between $`\eta _4`$ and $`\eta _2`$, namely

$$\eta _4=\frac{30\eta _2^211\eta _2+1}{16\eta _22}.$$
(24)
In appendix A, it is shown that this relation is precisely obeyed by the class of $`q`$– distributions Eq.(17) for which $`\stackrel{~}{\psi }(s)`$ was already solved.

The fact that the expansion of $`\stackrel{~}{\psi }(s)=[\stackrel{~}{p}^{}(s)]^{1/z}`$ only fails at fourth order implies that a mean field approximation, in which one explicitly assumes a product state, does give the exact results up to the third moment of $`p^{}(f)`$. This is precisely the reason why the mean field solution $`p^{\mathrm{mf}}(f)`$ differs only marginally from the real solution. To be more precise, the deviation $`p^{\mathrm{mf}}(f)p^{}(f)`$ should change sign $`4`$ times, since it does not affect all moments lower than $`f^4`$. A careful inspection of the numerical results in cop  for a $`q`$–distribution in which $`q=0.1`$ or $`q=0.9`$ shows that these small “wiggles” are indeed present. To magnify this effect, we show our simulation data in Fig. 2.

For $`z=3`$, the problems already appear at third order. Since we have $`\{3,0,0\}`$, $`\{2,1,0\}`$ and $`\{1,1,1\}`$, we encounter two independent equations for $`\psi _3`$. Again, it turns out that the equations can be solved if there is an additional relation between the $`q`$–moments:

$$\eta _3=\frac{15\eta _2^2\eta _2}{9\eta _2+1}.$$
(25)
For $`z>3`$, there are two independent third order equations as well, originating from $`\{3,0,0,0,\mathrm{}\}`$, $`\{2,1,0,0,\mathrm{}\}`$ and $`\{1,1,1,0,\mathrm{}\}`$. This problem can always be overcome by assuming a particular relation between the moments $`\eta _3`$ and $`\eta _2`$, corresponding to the special $`q`$–distributions of Eq.(17). Since at higher orders the number of equations per coefficient $`\psi _n`$ becomes increasingly high, there will be no other $`q`$–distributions than those of Eq.(17) that obey Eq.(14), and thus have an uncorrelated force distribution.

## IV Evolution of moments

Now that we know that, in general, correlations do exist in the stationary force distributions, it is interesting to study the nature of these correlations. In this section, we write the evolution of the moments as master equations, along the lines of Ref. jac . With this formalism, we will, in the next section, analyze the correlations by finding the stationary states of these master equations.

First, let us define the second moments of a distribution as

$$M_2(k)=f_if_{i+k}=𝑑\stackrel{}{f}f_if_{i+k}P(\stackrel{}{f}).$$
(26)
We have reintroduced the site–index $`i`$, and $`k`$ is a displacement–vector in a layer. As the system is translationally invariant, these second moments depend only on the displacement $`k`$. The recursion for these moments is obtained by combining Eq.(2) and Eq.(4) as

$`M_2^{}(k)`$ $`=`$ $`{\displaystyle \underset{\alpha ,\alpha ^{}}{}}\left({\displaystyle H(\stackrel{}{q})𝑑\stackrel{}{q}q_{j,\alpha }q_{j+k+\alpha \alpha ^{},\alpha ^{}}}\right)`$ (27)
$`\times M_2(k+\alpha \alpha ^{}).`$
Using the overline notation for the $`q`$–averages again, Eq.(27) becomes

$$M_2^{}(k)=\underset{\alpha ,\alpha ^{}}{}\overline{q_{j,\alpha }q_{j+k+\alpha \alpha ^{},\alpha ^{}}}M_2(k+\alpha \alpha ^{}).$$
(28)
This relation is reveals from which points (in correlation–space) the moment $`M_2^{}(k)`$ receives a contribution during a recursion step. However, it is in fact easier to consider the opposite relation, revealing how much a moment contributes to correlation space points during recursion. The “inverse” of Eq.(28) becomes

$$M_2(k)\overline{q_{i,\alpha }q_{i+k,\alpha ^{}}}M_2^{}(k+\alpha ^{}\alpha )\mathrm{for}\mathrm{all}\alpha ,\alpha ^{}.$$
(29)
This latter relation allows for a master equation type formulation, as we may write it in the form

$`M_2^{}(k)M_2(k)=`$
$`{\displaystyle \underset{\gamma }{}}W_\gamma (k\gamma )M_2(k\gamma )W_\gamma (k)M_2(k).`$ (30)
The transition rates are defined as

$$W_\gamma (k)=\overline{q_{i,\alpha }q_{i+k,\alpha ^{}}},$$
(31)
with $`\gamma `$ determined by the set $`\alpha ,\alpha ^{}`$ as

$$\gamma =\alpha ^{}\alpha .$$
(32)
In the current problem, where we consider second order moments, the transition rates are particularly simple. If $`k0`$, the $`q`$–averages are independent, and will always give the value $`1/z^2`$ (this only holds for $`q`$–distributions of the type Eq.(3)). If $`k=0`$, one encounters second moments of $`\eta (\stackrel{}{q})`$, as in Eq.(21). This leads to the following transition rates:

$`k=0`$ $``$ $`W_0(0)=\eta _2W_{\gamma 0}(0)={\displaystyle \frac{1z\eta _2}{z(z1)}},`$
$`k0`$ $``$ $`W_\gamma (k)={\displaystyle \frac{1}{z^2}}.`$ (33)
So, the moments evolve in an anomalous diffusion process, with differing transition–rates at the origin. For a detailed discussion of the corresponding dynamics, see Ref. jac . Note that this diffusion takes place in a $`d1`$ dimensional space, as $`\alpha `$, and therefore also $`\gamma `$, is a displacement in a layer. In the remainder of this paper we use the bold notation $`𝜸`$ whenever the displacement is really a vector.

The advantage of this somewhat formal representation is that we can take it over to higher order moments without further ado. The generalization of the master equation for the $`n`$th order moments $`M_n(𝐫)`$ becomes:

$`M_n^{}(𝐫)M_n(𝐫)=`$
$`{\displaystyle \underset{\gamma }{}}W_\gamma (𝐫𝜸)M_n(𝐫𝜸)W_\gamma (𝐫)M_n(𝐫),`$ (34)
with the position indices $`𝐫=(k_1,k_2,\mathrm{},k_{n1})`$, and the displacements $`𝜸`$ defined as

$$𝜸=(\alpha _1\alpha ,\alpha _2\alpha ,\mathrm{},\alpha _{n1}\alpha ).$$
(35)
The dimensionality of the diffusion process has now become $`(n1)(d1)`$. The transition rates can be calculated as

$$W_\gamma (𝐫)=\overline{q_{i,\alpha }q_{i+k_1,\alpha _1}\mathrm{}q_{i+k_{n1},\alpha _{n1}}}.$$
(36)
Analogous to the second moments, these transition rates are all $`1/z^n`$, as long as the indices of the position vector $`𝐫`$ are not equal to zero nor coincide. However, the differing rates make the problem complicated, because one has to deal with different transition rates at special points, lines, planes etc. in the space of diffusion.

One can now study the correlations at infinite depth by finding stationary states of the master equation for the moments. As a first attempt to construct a stationary solution, i.e. $`M_n^{}(𝐫)M_n(𝐫)=0`$, one can try a detailed balance solution. Detailed balance means that there is no flow of “probability” from one point to another. In that case, all terms of the sum on the right hand side of Eq.(IV) vanish individually, i.e.

$$W_\gamma (𝐫)M_n(𝐫)=W_\gamma (𝐫𝜸)M_n(𝐫𝜸)\mathrm{all}𝐫,𝜸.$$
(37)
This condition can also be formulated in terms of elementary loops, which are the smallest possible pathways from a point to itself. For all lattices in this study, these elementary loops are triangles, and we denote the three jump rates as $`(a,b,c)`$ or $`(a^{},b^{},c^{})`$ depending on the direction in which the loop is traversed. It is easily verified that the property

$$abc=a^{}b^{}c^{}$$
(38)
must be obeyed in all elementary loops in order to have a detailed balance solution. In the next section we show that correlations appear whenever the detailed balance conditions are not obeyed.

## V Higher order correlations

In this section, we study the nature of the correlations for $`q`$–distributions of the type Eq.(3) that do not fall into the special class of Eq.(17). We first solve the stationary master equation for the second order moments, for which we already know that there are no correlations (section III). For the triangular packing ($`z=2`$), correlations only show up at fourth order, and these fall off as $`1/r^5`$. For $`z3`$, there are third order correlations that also decay with a power law; for the fcc packing ($`z=3`$) the decay is $`1/r^4`$. Finally, we provide a simple relation to calculate the various exponents.

### V.1 Second order moments: no correlations

In order to get familiar with the structure of the master equations, we first consider the second order moments desribed by Eq.(IV). Away from the origin $`k=0`$, all transition rates of Eq.(IV) are identical. Therefore, the detailed balance condition Eq.(37) requires all $`M_2(k0)`$ to be identical. The value at the origin $`M_2(0)`$ has to obey a detailed balance condition for each $`\gamma 0`$, but these equations are identical for all $`\gamma `$ because the corresponding rates are the same. Putting $`M_2(k0)=1`$, one obtains the following stationary solution:

$$f_if_{i+k}=\{\begin{array}{cc}\frac{z1}{z(1z\eta _2)}\hfill & \text{}k=0\hfill \\ 1\hfill & \text{}k0\hfill \end{array}.$$
(39)
This solution precisely describes an asymptotic state without any 2–point correlations, as the average of the product $`f_if_j`$ equals the product of the averages for all $`ij`$. Of course, any multiple of Eq.(39), also forms a stationary solution of Eq.(IV). However, these solutions are physically irrelevant in the thermodynamic limit, where the lattice size $`\mathrm{}`$ jac . Moreover, we find that the asymptotic second force moment is solely determined by $`z`$ and $`\eta _2`$. For critical $`q`$–distributions one has $`\eta _2=1/z`$, leading to a diverging second moment.

### V.2 Third order moments

The diffusion of third order moments $`f_if_{i+k}f_{i+l}`$ takes place on a $`2(d1)`$–dimensional lattice, since there are two free parameters $`k`$ and $`l`$ of dimension $`d1`$. On this lattice, there are three special subspaces, namely $`k=0`$, $`l=0`$, and $`k=l`$, for which the transition rates of Eq.(36) differ from the bulk–value $`1/z^3`$. Moreover, the rates at the origin $`k=l=0`$ differ from both the bulk–rates and the rates on the special subspaces.

Let us first consider the triangular packing ($`z=2`$), for which the third order moments diffuse on a $`2`$–dimensional lattice, with differing rates on three special lines. As these lines are all equivalent, it is natural to draw them at an angle of $`120^{}`$, see Fig. 3. We then obtain a triangular lattice, with transitions to six nearest neighbors and two self jumps, which are “transitions” to the same lattice site ($`𝜸=\mathrm{𝟎}`$). The detailed balance condition between a special line and the bulk is naturally identical to the second order condition, implying the same ratio as in Eq.(39). As the transition rates at the origin are again identical for each $`𝜸\mathrm{𝟎}`$ (because of symmetry), one can construct the following detailed balance solution:

$$f_if_{i+k}f_{i+l}=\{\begin{array}{cc}\frac{\eta _2}{(12\eta _2)^2}\hfill & \text{, origin}\hfill \\ \frac{1}{2(12\eta _2)}\hfill & \text{, lines}\hfill \\ 1\hfill & \text{, bulk}\hfill \end{array}.$$
(40)
This means that there are also no $`3`$–point correlations for $`z=2`$: at the origin we encounter $`f^3`$, on the lines we have $`f_i^2f_{i+k}=f^2f`$, and in the bulk $`f_if_{i+k}f_{i+l}=f^3`$. It is easily checked that condition Eq.(38) is indeed satisfied in every elementary loop.

For the fcc packing ($`z=3`$), the third order moments diffuse on a $`4`$–dimensional lattice. Unlike the $`z=2`$ packing, it is not possible to construct a detailed balance solution in this case. First, we write the displacement vectors as $`𝜸=(\alpha ^{}\alpha ,\alpha ^{\prime \prime }\alpha )=(\gamma _1,\gamma _2)`$, where the $`\alpha `$’s and $`\gamma `$’s are 2–dimensional vectors (Fig. 1). One can jump away from the origin with two different rates, namely $`\overline{q_1^2q_2}`$ and $`\overline{q_1q_2q_3}`$. These rates correspond to $`\gamma _1=\gamma _2`$ (towards a special plane) and $`\gamma _1\gamma _2`$ (into the bulk) respectively. Checking the detailed balance condition in the elementary triangle origin–plane–bulk–origin, it turns out that Eq.(38) is only obeyed if $`\eta _3`$ and $`\eta _2`$ are related as in Eq.(25). Of course, this is precisely the case for the class of Eq.(17) for which we know that asymptotic factorization occurs. In general, however, it is not possible to construct a detailed balance solution for the third order moments. In the next paragraph, we show that the absence of detailed balance indicates that there are force–correlations, that decay with a power law; in this case the decay is $`1/r^4`$.

### V.3 Fourth order moments

The fourth order moments $`f_if_{i+k}f_{i+l}f_{i+m}`$ of the triangular packing diffuse on the bcc lattice depicted in Fig. 4. The three directions $`k,l,m`$ precisely define a bcc primitive cell ash . There are now differing rates at the origin as well as on lines and planes for which one or more indices coincide or are equal to zero. The precise geometrical structure is explained in appendix B. There are now two a priori different directions away from the origin, that is to corners $`f_i^3f_{i+1}`$ and to body centers $`f_i^2f_{i+1}^2`$. Checking the loop condition Eq.(38) for the loop origin–corner–body center–origin, one finds that it is only satisfied when $`\eta _4`$ and $`\eta _2`$ are related as in Eq.(24).

The question that emerges is: What are the stationary solutions of the master equation, when the detailed balance condition is frustrated at the origin? To answer this question we first consider a simplified version of the bcc–problem, as a first order approximation. In this simple version, we assume that all jump rates are $`1/z^4=1/16`$, except at the origin where we distinguish between the two different directions. Although we neglect the differing rates on the special lines and planes, the loop condition is still frustrated in the elementary loop origin–corner–body center–origin. Using $`_\gamma W_\gamma (𝐫)=1`$, we write the stationary master equation as

$$M(𝐫)=\underset{\gamma }{}W_\gamma (𝐫𝜸)M(𝐫𝜸),$$
(41)
or

$$\left[12W_0(𝐫)\right]M(𝐫)=\underset{\gamma 0}{}W_\gamma (𝐫𝜸)M(𝐫𝜸).$$
(42)
This allows us to eliminate the two self rates $`W_0`$ by means of a simple transformation:

$`\widehat{M}(𝐫)`$ $`=`$ $`\left[12W_0(𝐫)\right]M(𝐫),`$
$`\widehat{W}_\gamma (𝐫)`$ $`=`$ $`W_\gamma (𝐫)/\left[12W_0(𝐫)\right].`$ (43)
The sum over the new rates again adds up to unity and Eq.(42) becomes

$$\widehat{M}(𝐫)=\underset{\gamma 0}{}\widehat{W}_\gamma (𝐫𝜸)\widehat{M}(𝐫𝜸).$$
(44)
Hence we can omit the self jumps by first solving the equation for the “hatted” variables, and then transforming back to $`M(𝐫)`$. As $`M(𝐫)1`$ for large $`r`$, it is convenient to write

$$\widehat{M}(𝐫)=\frac{7}{8}\left[1+\delta \widehat{M}(𝐫)\right].$$
(45)
The quantity $`\delta \widehat{M}(𝐫)`$ is in fact the appropriate measure for correlations foot2 . After eliminating the two self rates, all jump rates have become $`1/14`$, except at the origin where the rates to the $`8`$ corners $`(c)`$ can differ from the rates to the $`6`$ body centers $`(b)`$. We therefore have

$$\widehat{W}_\gamma (𝐫)=1/14+\delta (𝐫)\epsilon _\gamma .$$
(46)
The rates to the corners are denoted by $`\epsilon _c`$ and those to the body centers by $`\epsilon _b`$. They fulfill the condition $`8\epsilon _c+6\epsilon _b=0`$. This results in the following equation:

$$\delta \widehat{M}(𝐫)\frac{1}{14}\underset{\gamma 0}{}\delta \widehat{M}(𝐫𝜸)=\frac{8}{7}\widehat{M}(\mathrm{𝟎})\underset{\gamma 0}{}\epsilon _\gamma \delta \left(𝐫𝜸\right).$$
(47)
Note that this is a discrete version of Poisson’s equation: the LHS is a discrete Laplacian and the RHS, originating from deviating rates, acts as a multipole around the origin. This equation is solved in appendix B by a Fourier transformation, leading to

$$\widehat{M}(𝐫)=\frac{7}{8}+\widehat{M}(\mathrm{𝟎})\underset{𝐤}{}\frac{E(𝐤)}{1D(𝐤)}\mathrm{exp}(i𝐤𝐫).$$
(48)
The functions $`D(𝐤)`$ and $`E(𝐤)`$ are defined in appendix B; $`1D(𝐤)`$ comes from the discrete Laplacian (in the continuum equation it would simply be $`k^2`$), $`E(𝐤)`$ is the Fourier transform of the source, and the sum over $`𝐤`$ is the inverse Fourier transformation running over the Brillouin Zone. The amplitude of the source $`\widehat{M}(\mathrm{𝟎})`$ can be obtained self–consistently, by setting $`𝐫=\mathrm{𝟎}`$. This involves a complicated integral over the Brillouin Zone (BZ) of the bcc–lattice; the outcome, however, will be of the order unity. The large $`r`$ behavior of the correlations is determined by the small $`𝐤`$ behavior, so $`E(𝐤)/(1D(𝐤))`$ has to be expanded around $`𝐤=\mathrm{𝟎}`$. The first term that gives a contribution is

$`{\displaystyle \frac{49\epsilon _c}{24}}{\displaystyle \frac{d𝐤}{V_{\mathrm{BZ}}}\frac{\left(k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2\right)\mathrm{exp}(i𝐤𝐫)}{k^2}}`$
$`{\displaystyle \frac{343\epsilon _c}{32\pi }}\left[5{\displaystyle \frac{x^2y^2+y^2z^2+z^2x^2}{r^9}}{\displaystyle \frac{1}{r^5}}\right].`$ (49)
The solution of Eq.(47) decays as $`1/r^5`$; the terms $`x^2y^2`$ etc. give the proper angular dependence. This result can be directly understood from the analogy with electrostatics. The solution of Poisson’s equation Eq.(47) can be expanded in asymptotically vanishing spherical harmonics: $`Y_{lm}/r^{l+1}`$. The symmetry of the bcc lattice allows only harmonics with $`l4`$, leading to the observed $`1/r^5`$ decay.

So we find that the stationary master equation for the moments becomes a discrete Poisson’s equation, and the presence of differing transition rates leads to a multipole source around the location of these rates, see Eq.(47). However, this source is only “active” if there is no detailed balance, since detailed balance leads to trivial solutions like Eq.(40) foot4 . Keeping this in mind, let us now investigate the real fourth order problem, including the differing rates at the special lines and planes. We argue that the asymptotic value is still approached as $`1/r^5`$, but the amplitude of this field will be modified. Since there is no detailed balance, the differing rates at the lines and planes will act as sources as well. Their amplitudes, however, will decay with increasing distance, since the “flow” associated with the absence of detailed balance becomes zero at $`r\mathrm{}`$. The effect of the induced sources at the special lines and planes can be taken into account perturbatively. The first step is to only consider the effect of the origin, as we have done above. The second step would be to compute the strength of the sources at the lines and planes on the basis of the first order solution, and then to determine their function $`E(𝐤)`$ and recalculate the solution Eq.(48). The induced sources around the origin basically lead to a modification of the strength $`\widehat{M}(\mathrm{𝟎})`$, but not of the asymptotic decay. However, the far away points at the lines and planes could modify the asymptotic decay. A closer inspection of the field of these sources shows that it is of order $`1/r^7`$, since the differing rates lead to a local Laplacian acting on the first order field decaying as $`1/r^5`$. Hence, every step of this perturbative calculation yields a leading term $`1/r^5`$; the amplitude changes in every step and its determination is a difficult problem indeed.

### V.4 Correlations for general $`z`$

From the previous section, it is clear that correlations occur whenever the detailed balance condition is frustrated around the origin. The stationary master equation then becomes a discrete Poisson’s equation in $`(n1)(d1)`$ dimensions, leading to correlations that decay with an integer power of the distance $`r`$. Following the derivation in appendix B, it is clear that the asymptotic behavior comes from the lowest non–isotropic term in $`E(𝐤)`$, since division by $`1D(𝐤)k^2`$ gives a singularity. The value of the exponent can be calculated as

$$(n1)(d1)+order2,$$
(50)
where $`(n1)(d1)`$ is the dimensionality of the correlation space and $`order`$ is the order of the lowest non–isotropic terms in the expansion of $`E(𝐤)`$. Although this result is remarkably simple, the actual calculation of $`E(𝐤)`$ is not trivial, as it reflects the symmetries of the jump directions on the $`(n1)(d1)`$ dimensional lattice. Working out the $`4`$–dimensional lattice of the third order moments in the fcc packing, we find that $`order=2`$ and correlations vanish as $`1/r^4`$.

## VI Correlated $`q`$–distributions

So far, we have only discussed $`q`$–distributions of the type Eq.(3), for which there are no correlations between $`q`$–values at different sites. We have shown that, for these $`q`$–distributions, there are no asymptotic $`2`$–point force correlations. In this section we will demonstrate that even the smallest correlation between $`q`$–values at different sites induces $`2`$–point force correlations. We first solve the problem for arbitrary correlations in the triangular packing. Then, we study the fcc packing assuming only a nearest–neighbor $`q`$–correlation; this already leads to force correlations that decay as $`1/r^6`$.

### VI.1 Triangular packing with arbitrary $`q`$–correlations

In general, the (second order) transition rates are defined by Eq.(31). For $`z=2`$, the displacement vector $`\alpha `$ can only take two values, for which we conveniently choose $`\pm \frac{1}{2}`$. This allows us to write the transition rates as

$`W_0(k)`$ $`=`$ $`\overline{q_{i,+\frac{1}{2}}q_{i+k,+\frac{1}{2}}}=\overline{q_{i,\frac{1}{2}}q_{i+k,\frac{1}{2}}},`$
$`W_{+1}(k)`$ $`=`$ $`\overline{q_{i,\frac{1}{2}}q_{i+k,+\frac{1}{2}}}`$
$`=`$ $`\overline{q_{i,\frac{1}{2}}(1q_{i+k,\frac{1}{2}})}=1/2W_0(k),`$
$`W_1(k)`$ $`=`$ $`\overline{q_{i,+\frac{1}{2}}q_{i+k,\frac{1}{2}}}`$ (51)
$`=`$ $`\overline{q_{i,+\frac{1}{2}}(1q_{i+k,+\frac{1}{2}})}=1/2W_0(k).`$
Asymptotically $`W_0(k)`$ has to approach the value 1/4, for $`q`$–distributions without long–ranged correlations. As the second moments diffuse on a line, one can easily construct a detailed balance solution:

$$\left[1/2W_0(k1)\right]M(k1)=\left[1/2W_0(k)\right]M(k),$$
(52)
or

$$M(k)=\frac{1/2W_0(0)}{1/2W_0(k)}M(0).$$
(53)
This is the general form of the $`2`$–point force correlations $`M(k)`$ in the triangular packing, as a function of $`W_0(k)`$ that describes the $`q`$–correlations. One can draw two interesting conclusions from this result. First of all, there can only be an uncorrelated solution if $`W_0(k)`$ is constant (i.e. $`1/4`$) for each $`k0`$. This means that even the smallest $`q`$–correlations lead to force correlations. Secondly, the long distance behavior of the $`2`$–point force correlations is identical to that of the $`2`$–point $`q`$–correlations, following from the simplicity of Eq.(53).

### VI.2 fcc packing with nearest neighbor $`q`$–correlations

Unfortunately, the analysis is much more complicated for the fcc packing, whose second order moments live on the $`2`$–dimensional triangular lattice of Fig. 5. We therefore allow only correlations between $`q`$–values at neighboring sites. Remember that one can easily construct an uncorrelated solution $`M(𝐫)`$ for uncorrelated $`q`$–distributions, Eq.(39), since all detailed balance conditions at the origin are identical by symmetry. This still holds when there are nearest neighbor correlations. However, the detailed balance condition will now be frustrated on the ring of surrounding sites, as these are connected in four a priori different directions, see Fig. 5. In analogy to the problem discussed in the previous section, the stationary master equation for $`\delta \widehat{M}(𝐫)`$ transforms into

$$\delta \widehat{M}(𝐫)1/6\underset{\gamma \mathrm{𝟎}}{}\delta \widehat{M}(𝐫𝜸)=\rho (𝐫).$$
(54)
The “charge density” $`\rho (𝐫)`$ is only non-zero around the frustrated ring, see appendix C. Again, it is a discrete version of Poisson’s equation, but now in $`2`$–dimensions. The solution can therefore be expanded in cylindrical harmonics, $`\mathrm{exp}(in\varphi )/r^n`$, and the six–fold symmetry of the lattice requires $`n6`$. The problem is again solved rigorously by Fourier transformation of Eq.(54). In appendix C we show that

$$\delta \widehat{M}(𝐫)\frac{\mathrm{cos}(6\varphi )}{r^6},$$
(55)
which is in agreement with the simple electrostatic picture.

So, for the fcc packing, we find that even a nearest–neighbor $`q`$–correlation leads to $`2`$–point force correlations that decay with a power law. This algebraic decay is generic for $`z3`$ since any $`q`$–correlations lead to a master equation whose detailed balance relations cannot be solved around the origin.

## VII Discussion

We have studied force correlations in the $`q`$–model at infinite depth, for general $`q`$–distributions. The calculated correlation functions are rather unusual: for $`q`$–distributions of the type Eq.(3), correlations only show up at higher orders, and these correlations decay with a power of the distance. The only exceptions are the $`q`$–distributions given by Eq.(17), which do produce a factorized force distribution. The results for the triangular packing and the fcc packing are summarized in Table 1. As an example, consider two different sites $`i`$ and $`i+k`$ in a layer of the triangular packing. Since there are no correlations in the second and third order force moments, we find $`f_if_{i+k}=1`$ and $`f_i^2f_{i+k}=f^2`$, independent of the distance $`k`$. However, the moments $`f_i^3f_{i+k}`$ and $`f_i^2f_{i+k}^2`$ are correlated and approach their asympotitic value as $`1/k^5`$. The fact that one has to go to higher orders to observe force correlations is the reason why numerical simulations only marginally differ from the mean field solutions cop . The (single–site) mean field solutions $`p^{\mathrm{mf}}(f)`$ are correct up to the third order moments, for the triangular packing. This implies that $`p^{\mathrm{mf}}(f)`$ “wiggles” around the real solution $`p^{}(f)`$; the deviation $`p^{\mathrm{mf}}(f)p^{}(f)`$ changes its sign $`4`$ times (Fig. 2).

Packings that have more than three $`q`$–values per site ($`z3`$) already have third order correlations. Also this time correlations only decay algebraically; for the fcc packing we find $`1/r^4`$. This algebraic decay can be understood from an analogy with electrostatics. The force moments evolve according to a master equation, and the corresponding stationary state is described by a discrete version of Poisson’s equation. The “source” turns out to be a multipole around the origin, which is only active whenever the master equation has no simple detailed balance solution. The moments therefore approach their asymptotic (uncorrelated) values algebraically. The value of the exponent depends on the dimension of the correlation space $`(n1)(d1)`$, and on the symmetry of the multipole, see Eq.(50).

Although in general correlations do exist, there is a special class of $`q`$–distributions, given by Eq.(17), for which there are no force correlations at all. This has been demonstrated by means of condition (14), which has a nice physical interpretation. It can be shown that the function $`\psi (s)`$ is the Laplace transform of the distribution of interparticle forces that live on the bonds connecting the particles: $`v_{i,\alpha }=q_{i,\alpha }f_i`$. Although the $`q`$’s leaving a site are correlated (they have to add up to $`1`$), the corresponding $`v_{i,\alpha }`$ can become statistically independent. It is only when this miracle happens that the force distribution becomes a product state. Nevertheless, the $`q`$–distributions for which this is the case range from infinitely sharp ($`r\mathrm{}`$) to almost critical ($`r0`$).

Finally, we found that there will be 2–point force correlations whenever the $`q`$–values of different sites are correlated. Even with only nearest neighbor $`q`$–correlations, the fcc packing has force correlations that vanish as $`1/r^6`$. Again, the triangular packing is less sensitive for correlations; the nature of the force correlations is identical to that of the $`q`$–correlations, Eq.(53).

Acknowledgements The authors would like to thank Wim van Saarloos, Martin van Hecke and Martin Howard for stimulating discussions.

## Appendix A Moments of $`q`$–distributions

This appendix is about the moments of the $`q`$–distri-butions, defined by

$$\overline{q_1^{n_1}q_2^{n_2}\mathrm{}q_z^{n_z}}=\eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}q_1^{n_1}q_2^{n_2}\mathrm{}q_z^{n_z}.$$
(56)
These different moments are not independent because of the $`\delta `$–constraint. As the distributions are normalised, the zeroth order moments are unity; the first order moments are, of course, all $`1/z`$. All second order moments, for which $`_in_i=2`$, can be described by only one free parameter. Defining $`\eta _n`$ as

$$\eta _n=\eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}q_1^n,$$
(57)
one finds

$`{\displaystyle \underset{i=1}{\overset{z}{}}}\overline{q_1q_i}=\eta _2+(z1)\overline{q_1q_2}=`$
$`{\displaystyle \eta (\stackrel{}{q})\delta \left(1\underset{\alpha }{}q_\alpha \right)𝑑\stackrel{}{q}q_1\underset{i=1}{\overset{z}{}}q_i}=1/z,`$ (58)
hence

$$\overline{q_1q_2}=\frac{1}{(z1)}(1/z\eta _2).$$
(59)
From a similar argument, one can derive for the third order moments

$$\overline{q_1^3}=\eta _3\overline{q_1^2q_2}=\frac{1}{(z1)}(\eta _2\eta _3).$$
(60)
For $`z=2`$ there is even a relation between $`\eta _3`$ and $`\eta _2`$:

$$1=\underset{ijk}{}\overline{q_iq_jq_k}=2\eta _3+6\overline{q_1^2q_2}\eta _3=\frac{3}{2}\eta _21/4.$$
(61)
For $`z=3`$, there is an additional third moment, namely

$`1={\displaystyle \underset{ijk}{}}\overline{q_iq_jq_k}=3\eta _3+18\overline{q_1^2q_2}+6\overline{q_1q_2q_3}`$
$`\overline{q_1q_2q_3}={\displaystyle \frac{1}{6}}(19\eta _2+6\eta _3).`$ (62)
The extension to higher orders and higher $`z`$ is straightforward.

For the special class of $`\eta (\stackrel{}{q})`$ defined in Eq.(17), one can calculate the moments $`\eta _n`$ from a generalization of Eq.(III) zinn

$$\eta _n=\frac{\mathrm{\Gamma }(zr)\mathrm{\Gamma }(r+n)}{\mathrm{\Gamma }(r)\mathrm{\Gamma }(zr+n)}.$$
(63)
In order to show that Eq.(24) is indeed obeyed by the special class (with $`z=2`$), we first invert Eq.(63) for $`n=2`$

$$r=\frac{12\eta _2}{4\eta _21}.$$
(64)
From this one can calculate $`\eta _4`$ as a function of $`\eta _2`$, which precisely results in Eq.(24). A similar inversion for $`z=3`$ leads to

$$r=\frac{13\eta _2}{9\eta _21},$$
(65)
from which one derives Eq.(25).

## Appendix B The bcc lattice

In the triangular packing, the fourth order force moments $`f_if_{i+k}f_{i+l}f_{i+m}`$ diffuse on the bcc lattice of Fig. 4, with differing jump rates on special lines and planes. In this appendix, we list these rates explicitly and we solve the corresponding stationary master equation.

The jump rates can be calculated from

$$W_\gamma (k,l,m)=\overline{q_{i,\alpha }q_{i+k,\alpha ^{}}q_{i+l,\alpha ^{\prime \prime }}q_{i+k,\alpha ^{\prime \prime \prime }}},$$
(66)
with the $`z^4=16`$ jump directions

$$𝜸=(\alpha ^{}\alpha ,\alpha ^{\prime \prime }\alpha ,\alpha ^{\prime \prime \prime }\alpha ).$$
(67)
As $`\alpha `$ can take the values $`\pm \frac{1}{2}`$, there are two self rates for which all $`\alpha `$’s are the same. As a consequence, there are $`14`$ outgoing directions, namely $`\pm (1,0,0)`$, $`\pm (1,1,1)`$, and $`\pm (1,1,0)`$ plus their permutations. The first two are directions for which three of the four $`\alpha `$’s are equal, and they correspond to the corners of Fig. 4; the third represents the jumps towards the body centers. If all position indices in Eq.(66) are different, the transition rates are simply $`1/z^4=1/16`$. On the special lines and planes where one or more position indices coincide, we encounter differing rates. The geometry of the problem is summerized in Table 2.

From this table we deduce the rates $`\epsilon _c`$ to the corners and $`\epsilon _b`$ to the body centers, which occur in relation (46). We find

$$\epsilon _c=\frac{\overline{q_1^3q_2}}{12\overline{q_1^4}}\frac{1}{14},\epsilon _b=\frac{\overline{q_1^2q_2^2}}{12\overline{q_1^4}}\frac{1}{14},$$
(68)
and one easily verifies from the property $`q_1+q_2=1`$ that the relation $`8\epsilon _c+6\epsilon _b=0`$ holds. In general, the rates do not obey the detailed balance condition Eq.(38) in the elementary loop origin–corner–body center–origin. Keeping only the rates in this loop as deviations from the bulk leads to equation (47). For the definition of the two functions $`E(𝐤)`$ and $`D𝐤)`$ we introduce two auxiliary functions: one for the contribution of the corners

$`\stackrel{~}{E}_c(𝐤)`$ $`=`$ $`{\displaystyle \frac{1}{4}}[\mathrm{cos}{\displaystyle \frac{k_x+k_y+k_z}{2}}+\mathrm{cos}{\displaystyle \frac{k_xk_y+k_z}{2}}`$
$`+\mathrm{cos}{\displaystyle \frac{k_x+k_yk_z}{2}}+\mathrm{cos}{\displaystyle \frac{k_xk_yk_z}{2}}]`$
and one related to the body centers

$$\stackrel{~}{E}_b(𝐤)=\frac{1}{3}(\mathrm{cos}k_x+\mathrm{cos}k_y+\mathrm{cos}k_z).$$
(70)
The two functions $`\stackrel{~}{D}(𝐤)`$ and $`\stackrel{~}{E}(𝐤)`$ are then given as

$`\stackrel{~}{D}(𝐤)`$ $`=`$ $`{\displaystyle \frac{4}{7}}\stackrel{~}{E}_c(𝐤)+{\displaystyle \frac{3}{7}}\stackrel{~}{E}_b(𝐤)`$
$`\stackrel{~}{E}(𝐤)`$ $`=`$ $`\epsilon [\stackrel{~}{E}_c(𝐤)\stackrel{~}{E}_b(𝐤)],`$ (71)
with $`\epsilon =8\epsilon _c=6\epsilon _b`$.

For the large $`𝐫`$ behavior we need the expansions for small $`k`$. One finds

$$\stackrel{~}{E}_c(𝐤)=1\frac{1}{8}k^2+\frac{1}{384}[k^4+4(k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2)]+\mathrm{}$$
(72)
and

$$\stackrel{~}{E}_b(𝐤)=1\frac{1}{6}k^2+\frac{1}{72}[k^4(k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2)]+\mathrm{}$$
(73)
From these expressions one derives the expansion

$`{\displaystyle \frac{\stackrel{~}{E}(𝐤)}{1\stackrel{~}{D}(𝐤)}}`$ $`=`$ $`{\displaystyle \frac{7ϵ}{24}}(1{\displaystyle \frac{7}{32}}k^2`$ (74)
$`+{\displaystyle \frac{7}{8}}{\displaystyle \frac{k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2}{k^2}}+\mathrm{}).`$
The first two terms in the expansion are regular and thus give rise to short range contributions. The last term leads to the asymptotic behavior, by means of the inverse Fourier transform

$$\frac{d𝐤}{V_{\mathrm{BZ}}}\frac{\left(k_x^2k_y^2+k_y^2k_z^2+k_z^2k_x^2\right)\mathrm{exp}(i𝐤𝐫)}{k^2}.$$
(75)
This integral can be evaluated by differentiation of the well–known

$$\frac{d𝐤}{V_{\mathrm{BZ}}}\frac{\mathrm{exp}(i𝐤𝐫)}{k^2}\frac{1}{4\pi r},$$
(76)
where a factor $`k_x`$ in Eq.(75) corresponds to applying $`/x`$. This leads to expression (V.3).

## Appendix C $`q`$–correlations in the fcc packing

In Eq. (54) we formulated the problem for the second moments in the fcc packing with nearest–neighbor $`q`$-correlations. The “charge density” $`\rho (𝐫)`$ on the right hand side of the equation is the product of the moment $`\widehat{M}(𝜸)`$, referring to the neighbors of the origin (all are the same by symmetry), with a function whose Fourier transform is given by

$$\stackrel{~}{E}(𝐤)=\underset{\gamma ^{},\gamma }{}w_{\gamma \gamma ^{}}\mathrm{exp}[i𝐤(𝜸+𝜸^{})].$$
(77)
The $`w_{\gamma \gamma ^{}}`$ are the deviations from the bulk transition rates $`1/6`$. These are only non–zero for the ring of nearest neighbors around the origin shown in Fig. 5:

$`w_0`$ $`=`$ $`\epsilon _0,w_1=w_5=\epsilon _1,`$
$`w_2`$ $`=`$ $`w_4=\epsilon _2,w_3=\epsilon _0+2\epsilon _1+2\epsilon _2.`$ (78)
The equalities reflect the symmetry of the triangular lattice. Inserting Eq.(77) into the Fourier transform of Eq.(54) leads to

$$\widehat{M}(𝐫)=2/3+\widehat{M}(𝜸)\underset{𝐤}{}\frac{\stackrel{~}{E}(𝐤)}{1\stackrel{~}{D}(𝐤)}\mathrm{exp}(i𝐤𝐫).$$
(79)
The consistency equation for $`\widehat{M}(𝜸)`$ follows by taking $`𝐫`$ as one of the nearest neighbors of the origin. The function $`\stackrel{~}{D}(𝐤)`$ is given by

$$\stackrel{~}{D}(𝐤)=\frac{1}{3}\left(\mathrm{cos}k_x+\mathrm{cos}\frac{k_x+\sqrt{3}k_y}{2}+\mathrm{cos}\frac{k_x\sqrt{3}k_y}{2}\right),$$
(80)
and $`\stackrel{~}{E}(𝐤)`$ can be expressed as

$$\stackrel{~}{E}(𝐤)/6=\epsilon _0[1\stackrel{~}{D}(2𝐤)]+2\epsilon _1[1\stackrel{~}{D}^{}(𝐤)]+2\epsilon _2[1\stackrel{~}{D}(𝐤)],$$
(81)
with the new function

$$\stackrel{~}{D}^{}(𝐤)=\stackrel{~}{D}(\sqrt{3}k_y,\sqrt{3}k_x).$$
(82)
For the asymptotic behavior of $`\widehat{M}(𝐫)`$ we must make an expansion of $`\stackrel{~}{E}(𝐤)/(1\stackrel{~}{D}(𝐤))`$. For the first two terms we find

$`{\displaystyle \frac{1\stackrel{~}{D}(2𝐤)}{1\stackrel{~}{D}(𝐤)}}`$ $`=`$ $`\mathrm{\hspace{0.17em}4}(1{\displaystyle \frac{3}{16}}k^2+{\displaystyle \frac{3}{256}}k^4`$ (83)
$`+`$ $`{\displaystyle \frac{1}{192}}{\displaystyle \frac{k_x^66k_x^4k_y^2+9k_x^2k_y^4}{k^2}}+\mathrm{}),`$
$`{\displaystyle \frac{1\stackrel{~}{D}^{}(𝐤)}{1\stackrel{~}{D}(𝐤)}}`$ $`=`$ $`3(1{\displaystyle \frac{1}{8}}k^2+{\displaystyle \frac{1}{128}}k^4`$
$`+{\displaystyle \frac{1}{288}}{\displaystyle \frac{k_x^66k_x^4k_y^2+9k_x^2k_y^4}{k^2}}+\mathrm{}),`$
and the third term is simply a constant. The asymptotic behavior is given by Fourier inversion of the first singular term in $`𝐤`$, i.e.

$`{\displaystyle \frac{d𝐤}{V_{BZ}}\frac{\left(k_x^66k_x^4k_y^2+9k_x^2k_y^4\right)\mathrm{exp}(i𝐤𝐫)}{k^2}}`$
$`{\displaystyle \frac{960}{\pi }}{\displaystyle \frac{\mathrm{cos}(6\varphi )}{r^6}}.`$ (85)
This integral can be obtained by differentiation of

$$\frac{d𝐤}{V_{BZ}}\frac{\mathrm{exp}(i𝐤𝐫)}{k^2}\frac{\mathrm{log}(L/r)}{2\pi },$$
(86)
where $`L`$ is the size of the system.


format: text

text:
# The Index Theorem for Homogeneous Differential Operators on Supermanifolds

## Prehistory

Israel Gelfand conjectured that since the index of any elliptic operator does not vary under deformations, it should be possible to express it in terms of the topological data of the problem. Atiyah and Singer were the first to prove this conjecture. Since then, there were given several proofs of this conjecture until finally Quillen not only noticed a manifestly “superish” nature of the index theorem — the index is the supertrace of an operator — but in two papers that drew a lot of attention to the then new field, supermanifold theory, gave an outline of the proof whose clarity was due to this nature. Getzler did a lot on elaboration and elucidation of the mechanism of this proof see \[Q\] and \[BGV\]. (Of particular interest are papers by Alvarez-Gomé and Witten, see \[GSW\], who link this, super, approach to physics.)

Quillen and his first followers, even Witten, only used the most superficial features of super-manifold theory: it was sufficient for the classical index theorem they were interested in. For the researchers of supermanifold theory proper it is tempting to attack a more general problem: to translate the index theorem and related topics, like the Riemann–Roch–Hirzebruch–Grothendieck theorem, from manifolds to supermanifolds. How to perform such a generalization is, however, quite unclear even after Manin’s contribution towards solution \[M\].

One of the obstacles: the de Rham cohomology of a (smooth) supermanifold coinsides with that of the underying manifold. The inadequacy of de Rham cohomology as an instrument for description of supermanifolds is manifest: the dimension of a vector superspace or a supermanifold is an element of the ring $`[\epsilon ]/(\epsilon ^21)=\epsilon `$, so attempts were made to label cohomology, as well as homology, by a pair of integers, rather than by a single one.

The two most successful approaches are those expressed in:

(1) papers by Manin and his students on homology of supergrassmannians: the idea to label Schubert supercells with elements of an analogue of the Weyl group. Several versions of the exposition was written; for a review see the translation \[M\] of VINITI’s collection. This review contains also an alternative and in some cases different description of super analog of the Weyl group, due to Penkov and Skornyakov. This rival approach, later improved yet further \[PS\], is especially good for superization of the Borel–Weil–Bott–Kostant-Bernstein–Gelfand–Gelfand theorem.

This approach, however, only embraced cell decomposition of quotients of simple Lie supergroups modulo parabolic subsupergroups, not the general supermanifolds. And nobody had yet succeded to introduce cohomology of supermanifolds labelled with two numbers except Shander:

(2) Shander’s integration theory concerning complexes of pseudoforms, see \[L4\], no. 33, Chapter 5.

Another obstacle for superization: one of the most conventional formulations of the index theorem — in terms of fixed points — baffles anyone who tries to formulate it on a supermanifold of purely odd dimension or for an odd differential operator. Again Shander is the only one who managed to see some ways out of this predictment, see \[SH\].

Summing up, one was tempted to conjecture that the index theorem on supermanifolds reduces to that on underlying manifolds. Indeed, such a theorem was soon published, see \[RS\].

I never believed in this theorem. The reason for my incredulity was based on two facts:

(a) still another formulation of a particular instance of the index theorem, due to Bott \[B\], in the case of an homogeneous, i.e., $`G`$-invariant, differential operator on the coset space $`G/H`$ states that the index theorem reduces to the Weyl character formula for the irreducible representations of the simple Lie group $`G`$;

(b) the superized Weyl character formula reduces to the representation of the underlying Lie algebra for typical representations only, otherwise it does not and the character of atypical representations is given by the so-called Bernstein–Leites character formula for $`𝔰𝔩(1|n)`$ \[BL\] and $`𝔬𝔰𝔭(2|2n)`$ \[L3\], respectively; for the most general form of the character formula see \[PS\], \[S\]. (Recall also earlier descriptions of the typical representation \[Be\] for $`𝔤𝔩`$ and $`𝔬𝔰𝔭`$; \[K\] for all simple finite dimensional Lie superalgebras with Cartan matrix; \[P\] for the series $`𝔮(n)`$ and the review \[L3\] for all representations, not only typical ones of the other series). The atypical representations are the most interesting ones: the identity one, the adjoint one, their symmetric powers, etc.

In \[L1\] I suggested to superize the index theorem along the lines of \[B\] but for various reasons the project was put aside. Recent discussions with B. Sternin and I. M. Gelfand encouraged me to reread \[B\] to see what can be sulvaged under superization; this note is a straightforward answer; conjecturally there exists an algebraic proof fitting any character formula.

I will confine myself to the case of unitary supergroups only. This, I am sure, is a purely technical restriction: the idea of the proof of the main theorem here follows word-for-word that of \[B\], where the compactness of $`G`$ is crucial. For the Lie groups this does not restrict the generality, since every simple Lie group has a compact form. Contrarywise, almost no simple complex Lie supergroup has a compact form and there are many compact superspaces $`G/H`$ with a noncompact $`G`$, see \[Se\]. Nevertheless, I am sure that the compactness of $`G`$ is beside the point here. The statement, being of purely algebraic nature, should be proved accordingly.

## Digest of Bott’s result

Let $`E`$ and $`F`$ be smooth complex vector bundles over a smooth compact manifold $`M`$; let $`D:\mathrm{\Gamma }(E)\mathrm{\Gamma }(F)`$ be an elliptic operator, see \[BGV\]. Set

$$\mathrm{ind}(D)=dim\mathrm{Ker}(D)dim\mathrm{Coker}(D).$$

It is now well-known, that if $`D`$ is elliptic, then $`\mathrm{ind}(D)`$ is finite (\[BGV\]).

Let $`M=G/H`$, where $`G`$ is a compact connected Lie group (in the proof we also need $`\pi _1(G)`$ to be without 2-torsion), $`H`$ a connected Lie subgroup of the same rank as $`G`$. Let $`E`$ and $`F`$ be induced by representations of $`H`$. We will say that $`D`$ is homogeneous if it is $`G`$-invariant. Passing from manifolds to supermanifolds we require the above for their respective underlying manifolds.

Define the ring $`R(G)`$ of virtual $`G`$-modules as follows. Let $`R(G)`$, as additive group, be generated over $``$ by the full set $`\{M_\lambda \}_{\lambda \mathrm{\Lambda }}`$ of irreducible $`G`$-modules. The representative of $`M_\lambda `$ in $`R(G)`$ is called its class and denoted by $`[M_\lambda ]`$. The ring structure of $`R(G)`$ comes from tensoring: $`[A][B]=[AB]`$. The map $`[]:M[M]`$ should possess the following properties of which only 4) is new as compared with \[B\]:

1) any additive function $`f:\text{G-mod}\text{Ab}`$ on the category of $`G`$-modules with values in the category of abelian groups Ab can be extended to a function $`f:R(G)\text{Ab}`$ such that $`f(M)=f([M])`$ and if the sequence $`0ABC0`$ is exact, then $`f([B])=f([A])f([C])`$;

2) any group homomorphism $`\phi :HG`$ extends to a ring homomorphism $`\phi ^{}:R(G)R(H)`$;

3) the map $`MdimM`$ (here $`dim`$ is the superdimension for modules over supergruops) can be lifted to a ring homomorphism $`dim:R(G)[\epsilon ]`$.

4) $`R(G)`$ is a $`[\epsilon ]`$-module with the action of $`\epsilon `$ given by the formula $`\epsilon [M]=[\mathrm{\Pi }(M)]`$.

Set $`\chi (D)=[\mathrm{Ker}D][\mathrm{Coker}D]`$. This is a refined index, since $`\mathrm{ind}(D)=dim\chi (D)`$.

Denote by $`\widehat{R}(G)`$ a completion of $`R(G)`$, the ring of formal power series on the same generators as $`R(G)`$. Let us extend the intertwining number of two $`G`$-modules

$$(A,B)=dim_{}\mathrm{Hom}_G(A,B)$$

to a symmetric bilinear form $`<,>_G`$ on $`R(G)`$ by setting $`<M_\lambda ,M_\nu >_G=\delta _{\lambda ,\nu }`$.

Given a group homomorphism $`i:HG`$, define the formal inducing $`i_{}:R(H)\widehat{R}(G)`$ by seting

$$i_{}:[A]\underset{\lambda }{}<i^{}(M_\lambda ,[A]>_H[_\lambda ],\text{ for any }AOb\text{H-mod},M_\lambda Ob\text{G-mod}.$$

For any homogeneous elliptic operator $`D:\mathrm{\Gamma }(E)\mathrm{\Gamma }(F)`$, where $`E`$ and $`F`$ are induced from $`H`$-modules $`K`$ and $`L`$, respectively, define the homogeneous symbol of $`D`$ to be $`s(D)=[K][L]R(H)`$.

###### Theorem .

1) $`\chi (D)=i_{}(s(D))`$.

1) For any $`G`$-module $`M`$ there exists a homogeneous elliptic operator $`D_M:\mathrm{\Gamma }(E)\mathrm{\Gamma }(F)`$ with $`\chi (D_M)=[M]`$.

The beauty of this theorem is twofold: the right hand side of 1) is the Weyl character formula of $`M`$ for $`D=D_M`$ and thanks to property 1) of the function $`[]`$ we see that $`[\mathrm{Ker}D][\mathrm{Coker}D]=[E][F]`$ does not depend on $`D`$! (Indeed, the sequence $`0\mathrm{Ker}DEF\mathrm{Coker}D0`$ is exact.)

To superize Bott’s result, we only have to replace in the above definitions groups with supergroups and dimension with superdimension, so the target ring in property 3) is $`[\epsilon ]/(\epsilon ^21)`$. To follow Bott’s proof, we need to integrate over a supermanifold with compact base which is only possible for the unitary supergroup. I am sure, however, that there is another, algebraic, proof of this theorem.


text:
# Large scale numerical investigation of excited states in poly(phenylene)

## Abstract

A density matrix renormalisation group scheme is developed, allowing for the first time essentially exact numerical solutions for the important excited states of a realistic semi-empirical model for oligo-phenylenes. By monitoring the evolution of the energies with chain length and comparing them to the experimental absorption peaks of oligomers and thin films, we assign the four characteristic absorption peaks of phenyl-based polymers. We also determine the position and nature of the nonlinear optical states in this model.

The phenyl-based conjugated polymers have attracted immense interest from physicists and chemists as a result of the discovery of electroluminescence in poly(phenylenevinylene) (PPV) . In particular, a great deal of theoretical and experimental effort has been devoted to understanding the important neutral excitations driving nonlinear optical processes, the positions of triplet states and the nature of exciton binding and decay . Correlation between $`\pi `$-electrons has proved to be important in obtaining an accurate description of excited states in conjugated systems. It is thus desirable to have accurate numerical solutions of correlated electron models for these systems (without recourse to uncontrolled approximation schemes such as configuration interaction approximation schemes, or perturbative schemes starting from a noninteracting molecular orbital or $`k`$-space band description). In this paper we provide, for the first time, accurate results for a realistic semi-empirical model for poly(phenylene) and its oligomers, a model system for phenyl-based conjugated polymers.

As our starting point for describing poly(phenylene), we adopt the Pariser-Parr-Pople (P-P-P) model ,

$``$ $`=`$ $`{\displaystyle \underset{<ij>\sigma }{}}t_{ij}[c_{i\sigma }^{}c_{j\sigma }+\mathrm{h}.\mathrm{c}.]`$ (3)
$`+U{\displaystyle \underset{i}{}}\left(n_i1/2\right)\left(n_i1/2\right)`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}V_{ij}(n_i1)(n_j1),`$
where $`<>`$ represents nearest neighbors, $`c_{i\sigma }`$ destroys a $`\pi `$-electron on conjugated carbon atom $`i`$, $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$ and $`n_i=n_i+n_i`$. We use the Ohno parameterisation for the Coulomb interaction, $`V_{ij}=U/\sqrt{1+(Ur_{ij}/14.397)^2}`$, where $`r_{ij}`$ is the inter-atomic distance in Å and $`U=10.06`$ eV . The transfer integrals ($`t_{ij}`$) and bond lengths are $`2.539`$ eV and $`1.4`$ Å respectively for phenyl bonds, and $`2.22`$ eV and $`1.51`$ Å respectively for single bonds . In the following we consider oligophenyls of $`N`$ repeat units (phenyl rings), as depicted in Fig. 3.

When solved accurately, the P-P-P theory has proved remarkably accurate in describing many of the low-lying excitations in a wide range of conjugated molecules . Given such debates as to the nature of exciton binding in PPV and the related issues of interpreting various spectral features , it is intriguing as to what the P-P-P theory predicts for large phenyl systems. Exact diagonalisation is currently restricted to systems with two phenyl rings. The advent of the density matrix renormalisation group (DMRG) method  has been fortuitous in that it has the potential to provide effectively exact numerical solutions of the P-P-P theory for systems with hundreds of conjugated carbon atoms. The DMRG has obstensively been restricted to systems with acetylene-like (simple, linear chain) geometry  to date. This is because it has not been possible to keep a complete basis for the phenyl repeat unit (4096 basis states) as well as sufficient system and environment block states. In this work we develop an accurate method for treating the phenyl block based on local Hilbert space optimisation .

Our approach starts with the assumption that an optimised local Hilbert space prepared in a small system should also be suitable for a large system. We employ the following procedure:

1. The P-P-P model is solved exactly for the 12-site ($`N=2`$) biphenyl molecule. The right-hand ring acts as the environment for the left-hand ring. A reduced density matrix for the left-hand ring is constructed by performing a partial trace of the target state(s) projection operator over the states of the right-hand ring. The $`m_A`$ optimal states for the left-hand ring (i.e. the optimised local Hilbert space) are those density matrix eigenstates with the highest eigenvalues. These are the natural many body states of the phenyl ring at the end of an oligo-phenylene molecule, in analogy to the natural one-electron orbitals of quantum chemistry.
2. For phenyl rings added to the middle of the lattice during the DMRG iterations (i.e. rings with neighbours to the left and right), we apply periodic boundary conditions to the biphenyl molecule to again construct an optimal basis, of which the $`m_B`$ most important states are retained. We note that this method of determining a suitable repeat unit Hilbert space from a small, exactly diagonalisable system is employed at the first stage of the “four-block” DMRG method .
3. Once the local Hilbert space has been optimised the usual inifinite lattice DMRG approach is used to build up oligophenyls. We use a three-block infinite lattice method, the environment block being the reflected system block. This enables us to apply the full symmetries of the Hamiltonian, which are the two planes of reflection, and particle-hole and spin-flip symmetries. We use the symbols $`A_g`$ and $`B_{1u}`$ to denote states that are symmetric under long-axis reflection, the $`A_g`$ states being symmetric under short axis reflection, the $`B_{1u}`$ states being antisymmetric under short axis reflection (polarised along the molecular axis). We use the symbol $`B_{2u}`$ to denote states that are polarised perpendicular to the molecular axis, i.e., states which are symmetric under long axis reflection and antisymmetric under short axis reflection. The spin-flip symmetry is used to distinguish singlet and triplet states, i.e., singlet states in the $`A_g`$ sector are denoted $`{}_{}{}^{1}A_{g}^{}`$ whilst triplet states are denoted $`{}_{}{}^{3}A_{g}^{}`$. Finally, the particle-hole symmetry of states is denoted by a “$`+`$” or “$``$” superscript. For example, the second lowest state in the $`A_g`$ singlet sector with positive particle-hole symmetry is denoted $`2^1A_g^+`$.
4. Oligophenyls with an even number of phenyl rings are constructed by performing DMRG iterations with no repeat unit (i.e., just with a system and environment block).

We embellish the basic approach described above with some further optimisations that serve to accelerate the convergence of energies with $`m_A`$ and $`m_B`$. Fig. 1 shows the construction of the initial phenyl bases: the end phenyl unit basis (for the first system block) is derived from the exact diagonalisation of biphenyl with open boundary conditions; the initial middle phenyl unit basis is derived from exactly diagonalising biphenyl with periodic boundary conditions. However, before using the middle block basis in DMRG iterations it pays to perform an improvement step using the trimer ($`N=3`$). Namely, we diagonalise the timer using a relatively small number $`m_A^{(0)}`$ of states in the system block, and a large number of states $`m_B^{(0)}`$ in the middle block. A reduced density matrix is then formed for the middle unit (by performing a partial trace over the system and environment block Hilbert spaces on an appropriate projection operator) and diagonalised to give the basis used for middle units in the subsequent DMRG iterations.

Fig. 2 shows a DMRG iterative step in the infinite lattice algorithm.

1. First, a system block and its reflection are combined without a middle unit in order to treat systems with an even number of repeat units. In this calculation $`m_A^{(1)}`$ states are retained in the blocks.
2. Next, a middle unit is inserted in the superblock. However, we perform two further basis optimisations to improve the middle block and system block bases in the context of this system size before obtaining our final energy estimates.
3. In the first step we improve the middle block basis (obtained from the trimer calculation depicted in Fig. 1) by again retaining a relatively small number $`m_A^{(2)}`$ states in the system and environment blocks and a relatively large number $`m_B^{(2)}`$ states in the middle block, with $`m_B^{(2)}m_B^{(0)}`$. A reduced density matrix is again formed for the middle block and diagonalised to give an improved middle block basis for this system size.
4. A similar optimisation is then carried out for the system block, this time keeping a relatively small number of states $`m_B^{(3)}`$ in the middle unit and a relatively large number of states $`m_A^{(3)}`$ in the system and environment blocks. A reduced density matrix is formed for the system block by tracing out over the middle unit and environment block spaces.
5. Finally, the energy estimates are obtained for this system size by retaining $`m_A<m_A^{(3)}`$ states from the reoptimised system block basis and $`m_B<m_B^{(2)}`$ states from the optimised middle block basis. It is also at this stage that the system block and middle unit are augmented to produce the system block used in the next iteration.

Before discussing our results, we assess the convergence of the DMRG procedure with $`m_A`$ and $`m_B`$. We first consider the noninteracting limit ($`U=0`$) where comparisons can be made with exact results. Table I shows the energies obtained for a number of states in various symmetry channels for $`N=11`$ phenyl rings using $`m_A=280`$ and $`m_B=468`$. In this example $`m_A^{(0)}=112`$ and $`m_B^{(0)}=606`$ and the extra middle unit and system block optimisations ((b) and (c) in Fig. 2) are not used. The accuracy varies with the state studied but even the most inaccurate state is resolved to within 0.02 eV. These results are a number of orders of magnitude more accurate than a previous effort to apply the DMRG to the oligo-phenylene geometry .

In the interacting case of interest the convergence of energies is monitored as a function of both $`m_A`$ and $`m_B`$. We first consider biphenyl ($`N=2`$) in order to get an idea of the efficiency of the phenyl bases. In Table II we list the energies of the two lowest $`{}_{}{}^{1}A_{g}^{+}`$ states as functions of $`m_A^{(1)}`$, the number of system block states retained in step (a) of the procedure depicted in Fig. 2. The indication is that even around 100 states around 3% of the total basis gives a reasonably accurate truncated basis for the phenyl block. We next consider the initial trimer step ($`N=3`$) depicted in Fig. 1. In Table III we list results for the $`2^1A_g^+`$ transition energy as a function of $`m_A^{(0)}`$ and $`m_B^{(0)}`$. The energy is converged to within around 0.01 eV. We have found that, once this initial reoptimisation of the middle phenyl block has been done for the trimer system, the parameter that controls the convergence (the parameter to which the energies are most sensitive) is the usual DMRG truncation parameter $`m_A`$, the number of states retained in the system and environment blocks. A typical example is given in Table IV where the $`2^1A_g^+`$ energy state is listed for various $`m_A`$ and $`m_B`$. The value of $`m_A^{(0)}`$ is set at 140 throughout, the value of $`m_B^{(0)}`$ is varied with $`m_B`$ as shown, and the extra middle unit and system block optimisations (steps (b) and (c) from Fig. 2) are not used. The results indicate that the $`2^1A_g^+`$ energy is resolved to within around 0.02 eV or better.

In production calculations we use all the embellishments depicted in Fig. 2. A typical example of the parameters used is given by: $`m_A^{(0)}=168`$, $`m_B^{(0)}=1000`$, $`m_A^{(1)}=340`$, $`m_A^{(2)}=140`$, $`m_B^{(2)}=740`$, $`m_A^{(3)}=290`$, $`m_B^{(3)}=300`$, $`m_A=180`$ and $`m_B=520`$. We have not carried out a comprehensive study of the effect of the embellishments (b) and (c) but, from running a few sets of parameters for the same target states, and from observing the change in energy obtained from using a reduced, reoptimised basis for either the system block or the middle phenyl unit, the indication is that they allow us to work with effectively larger values of $`m_A`$ and $`m_B`$. That is, using the embellished procedure with the above parameter set is probably equivalent to the simpler scheme (without steps (b) and (c)) with $`m_A290`$ and $`m_B740`$. We are fairly confident that, in the calculations in this paper, the lowest lying eigenstates in a given symmetry channel are resolved to within around 0.02 eV or better .

We now turn to a discussion of our results. We first embark on an investigation of the four characteristic absorption peaks of poly(phenylene) oligomers  and thin film polymers . Since the excitations of biphenyl ($`N=2`$) are relatively well understood, a knowledge of their evolution with chain length enables us to identify the key excitations in poly-phenyls. As the oligophenyls possess $`D_{2h}`$ symmetry, the dipole-active states are polarised either along the long axis ($`{}_{}{}^{1}B_{1u}^{}`$) or short axis ($`{}_{}{}^{1}B_{2u}^{}`$). Fig. 3 shows the $`N`$-dependence of transition energies of some key states, along with experimental results for biphenyl ($`N=2`$) (see ), oligomers ($`N=3,\mathrm{},6`$  and $`N=6`$ ) and polymer thin films ($`N=\mathrm{}`$) . We note that the particle-hole dipole-forbidden state ($`1^1B_{2u}^+`$) lies below the dipole-active $`1^1B_{1u}^{}`$ state in biphenyl. However, the $`1^1B_{2u}^+`$ state very weakly hybridises, and thus its energy is almost independent of chain length, converging to 4.4 eV. This energy agrees very well with the very weak second absorption peak at 4.4–4.5 eV  in polymer thin films. Adding weight to this interpretation is the observation of a weak but well defined 4.40 eV absorption peak in a highly textured film of sexiphenyl ($`N=6`$), orientated perpendicular to the substrate , as well as the weak, perpendicularly polarised absorption peak detected in orientated PFO film in the region 4.2–4.8 eV . In contrast, the $`1^1B_{1u}^{}`$ state strongly delocalises. As can be seen from Fig. 3, our DMRG results for the $`1^1B_{1u}^{}`$ in the $`N=3,\mathrm{},6`$ systems practically coincide with oligomer data . For large $`N`$ the $`1^1B_{1u}^{}`$ energy approaches 3.73 eV in reasonable agreement with the experimental peak observed at 3.63–3.68 eV , or 3.2 eV . Together these results indicate that the first (strong) and second (weak) absorption peaks in phenyl-based systems are the $`1^1B_{1u}^{}`$ and $`1^1B_{2u}^+`$ states, respectively.

The third and fourth absorption peaks are polarised normal ($`{}_{}{}^{1}B_{2u}^{}`$) and parallel ($`{}_{}{}^{1}B_{1u}^{}`$) to the long axis, respectively. In thin films they lie at 5.2–5.3 eV  and 5.7–6.0 eV  respectively. The $`{}_{}{}^{1}B_{1u}^{}`$ state is a localised intra-phenyl (Frenkel) excitation which lies at 6.16 eV in biphenyl . To track the position of this state as a function of $`N`$, we perform calculations targeting a number of $`{}_{}{}^{1}B_{1u}^{}`$ states together with the ground state $`1^1A_g^+`$ and calculate the dipole moments $`1^1A_g^+|\widehat{\mu }|j^1B_{1u}^{}`$, where $`\widehat{\mu }_ix_in_i`$ (with $`x_i`$ the $`x`$ coordinate of the $`i`$th atom) is the dipole operator along the long molecular axis. A typical set of dipole moments, for $`N=8`$, is given in Table V. In general, it is found that there are two $`{}_{}{}^{1}B_{1u}^{}`$ states that are strongly dipole connected to the ground state. The first is the $`1^1B_{1u}^{}`$, i.e., the first absorption peak. The second is the localised intra-phenyl excitation, i.e., the fourth absorption peak. (In the $`N=8`$ case this is the $`5^1B_{1u}^{}`$ state.) The $`N`$-dependence of this high-lying, strongly dipole connected $`{}_{}{}^{1}B_{1u}^{}`$ state is plotted in Fig. 3.

A conspicuous failure of the P-P-P model is its prediction for the remaining state, namely the lowest lying $`{}_{}{}^{1}B_{2u}^{}`$ state. The exact solution for biphenyl places this state at 6.66 eV, whereas experimentally it is at ca. 5.85 (and below the $`2^1B_{1u}^{}`$ state). As shown in Fig. 3 its calculated energy is 5.9 eV in the long chain limit, 0.6–0.7 eV higher than the experimental value. Recent work has shown that the predicted energy of this state is improved if dielectric screening is introduced into the Ohno interaction .

Fig. 3 also shows the lowest lying triplet ($`1^3B_{1u}^+`$) and the dipole-forbidden $`2^1A_g^+`$ state. Electroabsorption studies place the $`2^1A_g^+`$ state at around 4.6 eV , around 0.5 eV below the extrapolated P-P-P result of 5.1 eV. (A more recent study  places the $`2^1A_g^+`$ substantially lower). This discrepancy might possibly be explained in terms of the characteristic red shifts generally observed for certain excited states when going from well isolated chains to polymers in the solid state. Typical estimates for this polarisation or interchain screening shift are ca. 0.3 eV for the $`{}_{}{}^{1}B_{1u}^{}`$ state and ca. 0.6 eV for the $`2^1A_g^+`$ . As can be seen from Table V the $`2^1A_g^+`$ state has a large dipole moment with the $`1^1B_{1u}^{}`$ state, and unlike the case for polyenes, it is not pre-dominantly a pair of bound magnons, but a particle-hole excitation. (It is usually labelled the $`m^1A_g^+`$ state.) This particle-hole excitation is either a ‘p’-wave exciton, or the edge of the unbound particle-hole continuum.

To investigate this further, and the nonlinear optical properties of poly(phenylene) in general, we consider dipole moments between various $`{}_{}{}^{1}A_{g}^{+}`$ states and the $`1^1B_{1u}^{}`$ state, as well as between the $`2^1A_g^+`$ state and various $`{}_{}{}^{1}B_{1u}^{}`$ states. The $`N=8`$ values, listed in Table V, are representative of the general situation. We note that, in addition to the $`1^1A_g^+`$ and $`2^1A_g^+`$, another, higher lying state, which we denote the $`k^1A_g^+`$, also has an appreciable dipole moment with the $`1^1B_{1u}^{}`$. (For the $`N=8`$ case $`k=5`$.) We also observed a pattern in the $`2^1A_g^+|\widehat{\mu }|j^1B_{1u}^{}`$ values. Namely, the $`j=1`$ state has a strong dipole moment with the $`m^1A_g^+`$, as does the higher lying $`{}_{}{}^{1}B_{1u}^{}`$ absorption peak state (localised intra-phenyl exciton). In addition, there is another state, lying higher still, that has the largest dipole moment with the $`m^1A_g^+`$. We adopt the usual convention of denoting this state as the $`n^1B_{1u}^{}`$. (In the $`N=8`$ case $`n=7`$.)

In Fig. 4 we plot a number of $`{}_{}{}^{1}A_{g}^{+}`$ and $`{}_{}{}^{1}B_{1u}^{}`$ state transition energies as functions of $`1/N^2`$ . We see that there are a number of states in the $`{}_{}{}^{1}B_{1u}^{}`$ sector which seem to converge to the same energy as the $`1^1B_{1u}^{}`$ in the bulk limit. This is the band of ‘s’-wave charge-transfer excitons. Similarly, the $`m^1A_g^+`$ (i.e., the $`2^1A_g^+`$) has a number of momentum branches converging to the same energy, which is the band of ‘p’-wave charge-transfer excitons. Above this lies the localised intra-phenyl $`{}_{}{}^{1}B_{1u}^{}`$ state . Higher in energy still are the $`k^1A_g^+`$ and the $`n^1B_{1u}^{}`$ which appear to converge to the same energy in the $`N=\mathrm{}`$ limit. The strong dipole moments from the $`m^1A_g^+`$ to the $`n^1B_{1u}^{}`$ and the $`1^1B_{1u}^{}`$ to the $`k^1A_g^+`$, and the close proximity in energy of the $`k^1A_g^+`$ and the $`n^1B_{1u}^{}`$ states indicate that these mark the onset of the continuum of unbound particle-hole excitations. Lying below this continuum are the ‘s’- and ‘p’-wave charge-transfer excitons and the Frenkel exciton. This hypothesis could be checked further by calculating the average particle-hole separation for the various states to see if the $`n^1B_{1u}^{}`$ is the first unbound state in the $`{}_{}{}^{1}B_{1u}^{}`$ sector. The convergence of the $`k^1A_g^+`$ and the $`n^1B_{1u}^{}`$ energies to ca. 6.25 eV as $`N\mathrm{}`$ would imply a very large binding energy (ca. 2.5 eV) for the $`{}_{}{}^{1}B_{1u}^{}`$ exciton. However, band states are generally expected to be strongly effected by solid state screening (a red shift of 1.5 eV has been estimated for polyacetylene ). Taking such a shift into account would bring the $`n^1B_{1u}^{}`$ and hence the exciton binding energy much closer to the results implied by electroabsorption experiments .

In conclusion, we used a $`\pi `$-electron model to investigate the excited states of poly(phenylene), a model light emitting polymer. The transition energies were calculated using a new DMRG method. In order to accomodate the large repeat unit in the DMRG calculation, a truncated local Hilbert space was constructed for the phenyl unit. Using an optimal (or ‘natural’) basis affords us convergence in the excitation energies to within a small fraction of an eV. For the most part, the calculated transition energies agree well with the experimental results for biphenyl. By monitoring the evolution of the energies with chain length and comparing them to the experimental absorption peaks of oligomers and thin film polymers, we can assign the four characteristic absorption peaks of phenyl-based light emitting polymers. In ascending order they are the $`1^1B_{1u}^{}`$, $`1^1B_{2u}^+`$, $`1^1B_{2u}^{}`$ and a high lying, localised intra-phenyl (Frenkel exciton) $`1^1B_{1u}^{}`$ state. A failure of the P-P-P model is the incorrect prediction that the $`1^1B_{2u}^{}`$ state lies above the Frenkel exciton . We have also investigated the finite-size scaling of the nonlinear optical states in the longitudinally polarised one- ($`{}_{}{}^{1}B_{1u}^{}`$) and the two- ($`{}_{}{}^{1}A_{g}^{+}`$) photon sectors and found evidence for two $`{}_{}{}^{1}B_{1u}^{}`$ and one $`{}_{}{}^{1}A_{g}^{+}`$ excitonic bands lying below the continuum.

Computations were performed at the New South Wales Center for Parallel Computing, The Australian Center for Advanced Computing and Communications and The Australian Partnership for Advanced Computation. Dr Bursill was supported by The Australian Research Council and The J. G. Russell Foundation.