Dealing with long range interactions in the determination of polyelectrolyte ionization properties. Extension of the transfer matrix formalism to the full range of ionic strengths

The ionization state of charged macromolecules in solution is usually determined by the extent of the binding processes. These processes are very sensitive to the ionic strength of the medium, which are of long-range nature. The ionization properties of weak polyelectrolytes can be described by means of Ising-type models, which is only feasible when long-range interactions are neglected. Here, this formalism is extended to include long-range interactions by introducing a modified free energy involving only effective short-range interaction parameters. These parameters can be systematically calculated by using the GibbsBogoliubov variational principle. The technique is illustrated with the calculation of titration curves of homogeneous and heterogeneous polyelectrolytes in


Introduction
Binding of simple ions present in the background medium (protons, metal cations, etc.) controls the ionization state of weak polyelectrolytes in solution. This process is of paramount importance to understand the behaviour of charged macromolecules in a wide range of situations, from receptor-ligand interactions in biochemical systems [1,2], to the design of advanced coatings in materials science, wastewater treatment, pharmaceutical industry, etc. [3,4,5], and the role of natural organic matter in the geochemical cycling of trace metals [6], to mention just a few examples. Acid-base equilibria in weak polyelectrolytes represent the paradigmatic case of charge regulation in weak polyelectrolytes (due to the ubiquitous nature of proton ions in aqueous solution) and its study has been a classical topic in the past [7,8,9]. Among the different theoretical approximations used to describe the ionization properties in these systems, the site binding (SB) model has proven to be one of the most productive frameworks [10,11]. In this model, the ionization state of the molecule is defined as a set of protonating sites, which can adopt two possible states, i.e., protonated and deprotonated. For a linear polyelectrolyte, the resulting treatment is equivalent to the classical Ising model of ferromagnetism in the presence of an external field [12,13]. Potts models, a generalization of the Ising models which considers more than two possible states for the spins, have been applied to the calculation of conformational properties of linear polymers [14,15], coupling of ionization and conformational properties [16] or metal binding to polyelectrolytes [17].
The protonated sites experience both short and long range interactions. The latter are usually mediated by the solvent and can be described by simple potentials, such as the Debye-Hückel potential, or by suitable mean field approximations [18,19,20,21]. Short range interactions, however, can produce important correlations between the binding to neighbouring sites, so that the continuous mean field approach is no longer valid [22,23]. Moreover, they are mediated by the molecular skeleton rather than by the solvent, so that they cannot in general be modelled by simple interaction potentials. As a consequence, parameters accounting for chemically specific interactions must be included in the model [19,21,24,25]. Within the SB model, the problem can be tackled by expanding the free energy in a set of parameters (also known as "cluster" parameters) which involve interactions among an increasing number of sites: site protonation free energies, pair interaction energies, triplet (or three body) interaction energies, etc. [16,19,26,27,28,29,30,31]. In many cases, the resulting "cluster" expansion converges very fast to the exact free energy. If the number of sites is small (N ≤ 20), the necessary thermal averages can then be performed by direct enumeration. For systems with a large number of sites, direct enumeration becomes impracticable and other techniques, such as Monte Carlo (MC) simulation must be used [20,32,33,34,35,36,37,38,39,40,41,42].
For linear polyelectrolytes, the transfer matrix formalism is probably the most elegant and powerful method to compute thermal averages of Ising and Potts models. The key point is to relate the partition function of a system with N +1 units (spins, bonds or binding sites) to the partition of the system with N units. The resulting recurrence relationship can be expressed through the transfer matrix, whose elements are proportional to the conditional probability of finding the new unit in a particular state for a given state of the preceding one [12,43]. When applied to the binding of ions to polyelectrolytes, the method can be adapted to include a wide range of phenomena such as triplet interactions between sites [13], chelate complexation of metal ions [17], proton binding to polyampholytes [44,45], protein-DNA binding [46], super-capacitator charging [47], or coupling between ionization and conformational degrees of freedom [16,21,48]. This treatment, however, is only practicable if long range interactions are neglected, since the size of the transfer matrices grows exponentially with the range of the interactions [18]. In practice, this means that the study must be restricted to high ionic strengths, an important limitation specially for macromolecules which become unsoluble under such conditions [49,50,51,52].
In the present work, the transfer matrix method is modified in order to include long range interactions in an approximate but systematic way. With this aim, some necessary results on transfer matrices are recalled in section 2. The basic idea of the approach is presented in section 3. The original free energy is replaced by a new one involving only short range cluster parameters, which account for the long range interactions in an effective way. Equations for the new parameters are derived by using the Gibbs-Bogoliubov variational principle [43]. In the resulting formalism, both short and long range interactions are treated simultaneously. Although the results are valid for ionizable small molecules, surfaces or macromolecules of any architecture, they are particularly easy to implement for linear polyelectrolytes, since the transfer matrix method can be applied for the effective free energy in these systems as usual. The technique is illustrated with the calculation of titration curves of homogeneous and heterogeneous polyelectrolytes and polyampholytes in a wide range of ionic strengths (virtually, the method can be applied to the full range of ionic strengths achievable in practice). In the simplest form (which we will refer to as "first order correction"), the method computes effective site protonation free energies. The resulting isotherms reproduce with great accuracy the results obtained by MC simulations using the exact free energy. Corrections to other cluster parameters (higher order corrections), such as nearest neighbour interactions, are provided in section 4, although little significant improvement of the isotherms is obtained. The calculations are fast enough to perform straightforward fitting of binding parameters to experimental data. 3

Site binding (SB) model and transfer matrices
A particular protonation state (or microstate) of a system with N ionizable sites is characterized by a set of state variables s = {s i }where s i = 1 if the site is protonated and s i = 0 if the site is deprotonated. The microstate free energy can be expressed as the so called "cluster expansion" in the form pK i is the common logarithm of the protonation constant of the site i given all other sites are deprotonated, ij is the interaction energy of sites i and j, τ ijk represents the triplet interaction energy among sites i, j and k, and so on. The energy is expressed in thermal units, i.e., β = 1/k B T = 1, and the interaction parameters are divided by a factor ln 10 in order to be comparable in the pH scale. In most cases, the cluster expansion converges very fast to the exact free energy. Note that the choice of the protonation variables is not unique. For instance, for polyanions such as polycarboxylic acis, it could seem more natural to work with the variables q i = 1−s i , where q i = 1 where the site is charged (i.e. deprotonated) and q i = 0 is the site is uncharged (i.e. protonated). However, it has been shown that the choice of the protonation variables does not affect the resulting thermal averages if the cluster parameters are properly redefined [19]. The foregoing arguments and calculations are thus applicable regardless of the charge sign of the protonated sites. It is also worth mentioning that although the flexibility of the chain is not explicitly included in F (s), it can be shown that the cluster expansion (1) is still exact if the cluster parameters are regarded as proper averages over the conformational degrees of freedom [16]. The probability of a particular microstate in the semi-grandcanonical ensemble is given by where n is the total bound protons, a H is the proton activity, and Ξ = s e −F (s) a n H is the semi-grandcanonical partition function. For the purposes of this work, it is convenient to work in terms of the reduced free energy where µ i = (pH − pK i ) = log (K i a H ) is the reduced chemical potential in the pH scale. The relevant physical properties can be obtained by performing the proper thermal averages, which in many cases can be expressed in terms of derivatives of Ξ, related to the reduced free energy as For instance, the average degree of protonation of a particular site i is obtained by where Ω = − ln Ξ is the thermodynamic potential or free energy associated to the semi-grandcanonical ensemble. The average number of bound protons is given by As we will see later, it is also important to determine the correlation of the protonation degrees of two sites i and j, which can be expressed as The computation of the partition function can be performed by direct enumeration of the microstates only for N ≤ 20. Otherwise, methods borrowed from statistical mechanics must be used. For 2D systems, only a few cases can be exactly solved [12]. For 1D systems, however, transfer matrix techniques can be applied to a wide range of models, although only if long range interactions are neglected. The transfer matrix method consists of relating the partition function for a system with N+1 sites with that of a system with N sites in a recursive way. The link between both partition functions can be shown to be the transfer matrices. As a result, the partition function can be expressed as where T i is the transfer matrix of site i, whose elements T i ; σ,ρ are the Boltzmann factors corresponding to the increase in the reduced free energy. They are proportional to the conditional probability of adding a site i in a state ρ provided that the previous site i -1 was in the state σ. q and p are proper initiating and terminating vectors which depend on the details of the end of the chain. For identical sites eqn. (8) reduces to where T = T 1 = T 2 = · · · = T N . In the limit N → ∞, Ξ becomes proportional to λ N , the maximum eigenvalue of T , so that Ω ∼ −N ln λ and the average degree of protonation θ = θ i = ν/N reads A useful altervative expression to (10), which involves derivatives of the transfer matrix instead of λ, is derived in the appendix. The result is If Q is the matrix whose columns are the eigenvectors of T, a = (Q 1r , Q 2r , Q 3r , · · · ) represents the eigenvector corresponding to λ , and b is the r th row of An equivalent expression for the correlation function h ij is given by In the simplest case, when only nearest neighbour interactions are considered, the transfer matrix adopts the form [19] where z =Ka H is the reduced activity and = − log u is the interaction free energy between consecutive sites. The corresponding initial and final vectors are p = (1, 1) and q = (1, 0), respectively. If next-nearest neighbour and triplet interactions are included, the transfer matrix can be shown to be [45] T = where ζ ≡ i,i+2 = − log v represents the next-nearest neighbour interaction free energy and τ = − log w the triplet interaction free energy, only present when three adjacent sites are protonated. In this case, p = (1, 1, 1, 1) and q = (1, 0, 0, 0).

Effective site protonation free energies: first order correction
If long range interactions are included in the free energy, approximate methods are necessary. When the range of the interactions is increased in one site, the size of the transfer matrix is multiplied by a factor of 2. For instance, in order to include next-nearest neighbour approximation we need the 4 × 4 transfer matrix (14); 8 × 8 matrices are needed to account for next-next nearest neighbour interactions, and so forth [18]. Therefore, the size of the matrices grows exponentially with the range of the interactions and the method becomes impractical. The approach here proposed to overcome this difficulty is to replace the original free energy by an approximate one which only includes short range interactions, but which accounts for the long range interactions in an effective way. With this aim, the Gibbs-Bogoliubov variational principle [43] will be used. Let us assume that H (s) can be expressed as a sum of two terms so that H 0 (s, α) can be exactly solved, and α = {α k ; k = 1 · · · P } is a set of P parameters, not necessarily present in H (s). α k are chosen by means of some physical criterion such that ∆H is "small" compared to H 0 . Then, it can be shown that if H (s) is replaced by the new microstate reduced free energy H = H 0 + ∆H 0 , the Gibbs-Bogoliubov inequality holds where Ω 0 (α) = − ln Ξ 0 and · · · 0 represent the free energy and the thermal average corresponding to H 0 , respectively. Ω is thus an upper bound of Ω, so that the best approximation to Ω is obtained by choosing α k such that Ω is minimum The thermal averages are then approximated by using H (s, α opt ) instead of H (s), which in turn coincide with the thermal averages calculated using H 0 (s, α opt ), which is exactly solvable

Identical sites
Let us consider a linear chain composed by N identical but interacting sites. The protonated sites experience short range interactions with m neighbouring sites and long range interactions with the rest of sites by a potential φ ij . Let us also assume that there is no interaction between different chains. The reduced free energy is expressed as the sum of two terms The first brackets contain the short range interations and the second ones the long range interactions. Let us split H into two parts H = H 0 + ∆H defined as In doing so we try to find a correction µ to the site protonation reduced free energy, so that the main effects of the long range interactions are embedded in H 0 . µ can also be understood as a correction to the site pK -value, so that pK eff = pK − µ is the effective pK -value. The exact free energy Ω is then replaced by is the long range energy averaged over the unperturbed system, whose correlation function (7) is h 0 ij . Taking ∂ Ω/∂µ = 0, one obtains for µ where ν 0 = i s i 0 = (∂Ω 0 /∂µ ) / ln 10 = (∂Ω 0 /∂µ) / ln 10 is, according to (6), the average number of protons for the unperturbed model. An alternative to eqn. (22) is to determine µ by directly minimizing (??) using a suitable optimization routine. Eqn. (22), however, provides a transparent interpretation for µ , which results to be the average change of long range interaction energy involved in adding a new proton to the chain (calculated using H 0 ). Note that no consideration has been made about the shape or dimensionality of the system, so that eqn. (22) is valid for small molecules, surfaces or polyelectrolytes. For surfaces, not much is known about Ω 0 for arbitrary µ−values [12,53]. However, in problems involving linear chains, short range interactions can always be treated by means the suitable transfer matrix, with the reduced activity z replaced by where µ is calculated by using (21)(22). For N → ∞, ϕ becomes which can be evaluated by means of eqn. (12). If only nearest neighbour interactions are considered, the maximum eigenvalue λ 0 of the transfer matrix (13) can be analytically evaluated and, using (10) or (11), the average degree of protonation θ 0 = ν 0 /N remains [30] If the correlation between the long range interacting binding sites is neglected, we can assume that for j > i + m Introducing (26) in (22) we obtain which together with (23)(24)(25) generalizes the Frumkin isotherm, since it includes long range interactions within the mean field approximation, while short range interactions are treated in an explicit way. However, due to the drastic approximation (26), eqns. (25)(26)(27), unlike eqns. (22)(23)(24), are not accurate for the whole range of ionic strengths.     1a shows the titration curves corresponding to a linear chain with identical sites and nearest neighbour interactions at the short range level. As discussed above, short range interactions are mediated by the molecular skeleton rather than the solvent [19], so that specific values for them are necessary. Long range interactions are described by the Debye-Hückel potential for an extended linear chain where B 0.7 nm is the Bjerrum length in water at 298 K , b is the separation between consecutive protonating sites and κ −1 (nm) = 0.304/ I (M) is the Debye length. Note that the choice of the DH potential implies that the distance of two interacting sites is constant so that the chain is assumed to be rigid. Consequences of the chain flexibility in the ionization properties are analyzed, for instance, in [16,21,55] and in [56,57,58,59] The chosen parameters are pK = 9, b=0.2 nm and = 2. They are typical values for polyethylene(imine) [16,54]. Dots represent the values obtained from MC simulations, while continuous lines correspond to those obtained by using eqns. (21)(22)(23)(24)(25). The agreement between theory and MC simulations is excellent. Surprisingly they can be regarded as exact from the practical point of view. It is worth noting that the variation of the profiles reveals the importance of including the long range interactions even for relatively high ionic strengths (0.1 M), being completely necessary at low ionic strengths. In Fig. 1b the correction to the protonation free energy, µ is depicted. As expected, it increases in lowering the pH (increase in the charge), and in decreasing the ionic strength (lower electrostatic screening). Note that the wavy behaviour at pH values for which θ 0.5 corresponds to the maximum correlation between protonated and deprotonated sites, which tend to appear alternated in the chain [30]. As will be shown in section 4, this behaviour desappears when higher order corrections are included, since part of the effect is embedded in the effective nearest neighbour interactions. Fig. 2 shows the titration curves corresponding for the same model, but now triplet interactions and next-nearest interactions have been added at the short range level. The chosen parameters are τ = 0.43 while ζ has been calculated for each curve using the Debye-Hückel potential (28) taking |j − i| = 2, and introduced in the transfer matrix (14), so that correlations between next-nearest neighbour sites are treated in an explicit way. Triplet interactions are necessary to explain the titration curves of some small molecules and polyelectrolytes [13,19,28,30,54]. Recent papers have reported that they can be understood as the result of the coupling of conformational and ionization degrees of freedom [16,21]. Dots represent MC simulations, while continuous lines come from eqns. (21)(22), (23) and (10), together with the transfer matrix (14). A very good agreement is again obtained. Finally, we would like to comment that proper modifications of the free energy (19) could account for the interaction between two or several chains, by including the suitable interacting terms. Interaction between polyelectrolytes could drive interesting phenomena, such as charge regulation or phase separation, previously reported in colloidal suspen- sions [23,60,61,62,63]. This important topic is, however, out of the scope of this work.

Heterogeneous polyelectrolytes
The treatment presented in the previous section can be straightforwardly adapted to molecules with arbitrary shape and chemical composition, by correcting each site protonation free energy pK i , i = 1, · · · N , so that pK eff i = pK i − µ i . As an illustrative example, let us consider a linear polyelectrolyte with two different kind of sites, A and B, such that they form the alternating structure ABABAB... The corresponding reduced free energy is given by i=2,4,6,...
We need corrections to both site protonation free energies, µ A and µ B , and (30) in (16) the trial free energy Ω is given by  (31) with respect µ A and µ B and setting to zero, we obtain have been used. Again, we can numerically solve the system of equations (32) or, alternatively, find µ A and µ B by direct optimization of Ω using a suitable routine.
Eqn. (32) can be modified such that it provides a direct physical interpretation for µ A and µ B . Eqn. (33) states that, for given values of µ A and µ B , there corresponds a unique pair of values of ν 0 A and ν 0 B . Let us now suppose that we work in terms of the new variables ν 0 A and ν 0 B instead of µ A and µ B . Since the matrix J 1 is the jacobian of the transformation (33), eqns. (32) for µ A and µ B can be inverted and rewritten in terms of ν 0 A and ν 0 B in the simpler form which states that µ A (µ B ) can be interpreted as the change in the long range free energy ϕ when a new A (B) site is protonated, keeping constant the number of protonated B (A) sites. For linear polyelectrolytes H 0 can be solved using suitable transfer matrices. For the alternating sites model T is given, as in (13), by where z eff i = z i 10 −µ i are the effective reduced activities, corrected using eqns. (32)(33), and = − log u represents the nearest neighbour interaction energy between A and B sites. The resulting titration curves are compared to MC simulations in Fig. 3a. with parameters pK A = 7, pK B = 9, = 1.5, and b=0.2 nm, which represent alternating copolymers with two different basic functional groups. The agreement is very good for the full range of ionic strengths. In Fig. 3b the corrections to the site protonation energies are plotted as a function of pH. As expected, they increase in decreasing the ionic strength due to the lowering of the screening of coulombic forces. At ionic strengths lower than 0.1M, the correction in the protonation free energy at θ 0.5 becomes larger than 1 pK -unit and the effect of the long range interactions cannot be neglected. Note that the correction is larger for B sites than for A sites, since the latter present a lower affinity for the proton.
It is also interesting to test the method if, not only repulsive, but also attractive interactions are present. Attractive interactions are expected to produce high correlations between groups with opposite charge so that their presence is a good test to the present treatment. Let us consider a polyampholyte composed by two kind of alternating groups, forming a structure ABABAB... A groups are acidic (for instance, carboxylic groups) and negatively charged at high pH-values, while they are not charged at low pH-values. B groups have basic character (for instance, amino groups) and they are positively charged at low pH-values, while they are uncharged at high pH-values. The corresponding transfer matrix is where = − log (u) represents the interaction between a negatively charged (i.e. deprotonated) A group and a positively charged (i.e. protonated) B group, so that < 0 . Fig 4a compares the theoretical and simulated titration curves, for a polyampholyte with parameters pK A = 4.5 , pK B = 6, b= 0.2 nm and = −1.5. Again a very good agreement is obtained for the full range of ionic strengths, a particularly important point when working with polyampholytes since they can become insoluble at high ionic strengths [49]. It is worth to note that the isoelectric point, which for this particular model coincides with θ = 1/2, is independent of the ionic strength. Fig 4b presents the average protonation degree of sites A and B, θ A and θ B , together with the average charge per monomer (q) and the number of zwitterions (D = i s i s i+1 ) per monomer. The isoelectric point appears around pH 5, which, as expected, coincides with the maximum number of zwitterions. The effective parameters capture both the strong binding correlations due to the short range interactions and the more delocalized effect of long range interactions.      expected to converge very fast to the exact free energy. For simplicity, let us firstly treat the situation in which the short range interactions are described by the nearest neighbour model. Besides the effective reduced chemical potential, we also search an effective interaction energy. The free energy can now be splitted into The trial free energy adopts the form where D 0 = i s i s i+1 0 is the number of neighbouring protonated sites over the unperturbed system. Minimizing (38) with respect to µ and we obtain where has been used. Note that, as pointed out in subsection 3.2, Ω is the Legendre transformation of Ω 0 corresponding to the change of variables (40). A physical interpretation of µ 2 and arises using an argument similar to the one used in section 3.2. Eqn. (40) can be seen as the transformation between the two pair of variables (µ 2 , ) and (ν 0 , D 0 ) and J 2 is the jacobian of such a transformation. Inverting (39) , µ 2 and remain so that µ 2 represents the average increase in the long range energy when a new proton is bound (keeping constant the number of neighbouring interactions), while accounts for the change in the long range energy in creating a new neighbouring interaction (at a constant number of bound protons). This result explains why the first order correction works so well. Intuitively, one can expect that the change in the long range energy in creating a new neighbouring interaction is very small if no new proton is bound to the molecule. This is confirmed by the calculations, which show that the second order correction does not add significant improvement to the titration curves. In Fig. 1a (homogeneous polyelectrolyte) and Fig. 4a (polyampholyte) the second order correction (39) has been plotted (crossed markers). No significant improvement is obtained compared to the ones calculated by using the first order correction. In fig. 5 the second order corrections µ 2 and are plotted versus pH, for the same homogeneous polyelectrolyte depicted in Fig. 1, and compared to µ , the first order correction to the site protonation free energy. As expected, is small compared to µ 2 . We also observe that the wavy behaviour of µ is no longer present in µ 2 , but it has been replaced by the contribution to the nearest neighbour interaction, , which is maximum at pH-values for which θ 0.5. This could be related to the increment in free energy produced when the structure of alternating protonated and deprotonated sites is modified when two protonated sites become neighbours. Using the same procedure, higher order corrections to the triplet, next-nearest neighbour interactions can be calculated, and expressions of the type (41) derived. Similarly, the correction to a particular cluster parameter can interpreted as the average change in the long range energy when a new interaction (of the type described by the cluster parameter), is created in the polymer. However, as discussed above, no gain in accuracy is obtained. Therefore, at least for the models here studied, the first order correction to the protonation free energy is enough to reproduce almost exactly the MC simulations.

Conclusions
An approach to Ising-type models in polyelectrolytes which treats short and long range interactions simultaneously, is presented. The effect of long range interactions is included by means of effective short range cluster parameters, which can be systematically calculated by means of the Gibbs-Bogoliubov variational principle. The first order correction to the site protonation free energy correspond to the change of the average long range interaction energy when a new proton is added, a result which is valid for small molecules, polyelectrolytes and surfaces. However, it is specially straightforward to apply to linear molecules, where the transfer matrix method can be implemented for a wide variety of situations. In general this formalism has been used when long range interactions are neglected, since the size of the matrices grows exponentially with the range of the interactions. This fact often restricted the study of the titration curves to high ionic strengths, sometimes a hard limitation since many molecules of interest become unsoluble under such conditions. This difficulty is overcome with the treatment here presented.
Different models have been used to test the method. In all the cases investigated, where long range interactions have been included by means of the Debye-Hückel potential, the first order correction provides results which are almost indistinguishable from those obtained from MC simulations, allowing fast parameter fitting to experimental data for the full range of ionic strengths. As examples of heterogeneous systems, polyelectrolytes and polyampholytes with two kind of alternating sites have been chosen to test the method. In the case of polyampholytes, both repulsive and attractive interactions are present. Effective protonation free energies have been calculated for each chemically different site. It is important to emphasize that long range interactions imply important corrections to the titration curves even for not too low ionic strengths, which could lead to a significant bias in the fitted parameters if this type of interactions is not take into account.
Higher order corrections have also been implemented. In the second order correction, both the site protonation free energy and the nearest neighbour interaction are recalculated to include long range effects. It is shown that the correction to the nearest neighbour interaction energy corresponds to the average change in the long range energy when a new neighbouring proton-proton pair is created (keeping constant the number of bound protons). This justifies why the second order correction is so small and add little improvement to the resulting titration curves. As a general statement, we can say that the correction to a particular cluster parameter represents the average change in the long range energy when a new interaction (of the type described by the cluster parameter), is created in the polymer.
The basic ideas here presented could possibly be extended to a wider group of situations where Ising and Potts models are useful, such as ionization of surfaces and small molecules, conformational properties of macromolecules, polymerpolymer and polymer-surface interactions, molecular stretching, among others.
APPENDIX: degree of protonation and binding correlation in terms of derivatives of the transfer matrix for N → ∞ For a linear chain of N +1 identical sites the partition function is given by eqn. 9. In the limit N → ∞ the average degree of protonation θ is the same for all the binding sites, so that we can choose, for instance, the site located in the middle of the chain. θ is thus given by where N is even and T can be expressed in the diagonal form T = QΛQ −1 , where Q is a matrix whose columns are the eigenvectors of T, and Λ is a diagonal matrix. Let us suppose that the maximum eigenvalue is the (r, r) element of Λ and the corresponding eigenvector is the r th column of Q and Λ. where E is a matrix whose elements are zero excepting E rr = 1. Therefore, eqn. where a = (Q 1r , Q 2r , Q 3r , · · · ) contains the eigenvector corresponding to λ and b is the r th row of Q −1 , b = Q −1 r1 , Q −1 r2 , Q −1 r3 · · · . Introducing (A.6) in (A. 5) and taking into account that qa T and bp T are scalar quantities, the sought result (11) is obtained An expression for the correlation function h ij = s i s j can be derived in a similar way. For N → ∞, the correlation function h ij becomes independent of i, so that h ij can be expressed