text:
# 1 Introduction

## 1 Introduction

Let $`\tau _n`$ denote the maximal kissing number in dimension $`n`$, that is, the greatest number of $`n`$-dimensional spheres that can touch another sphere of the same size. Although asymptotic bounds on $`\tau _n`$ are known , little is known about explicit constructions, especially for $`n>32`$. Up to now the best explicit constructions have come from lattice packings. The kissing number $`\tau `$ of the Barnes-Wall lattice111The subscript gives the dimension. $`BW_n`$ in dimension $`n=2^m`$ is $`_{i=1}^m(2^i+2)`$, although for $`m5`$ this is weak (146,880, 9,694,080 and 1,260,230,400 in dimensions 32, 64 and 128, for example). In contrast, Quebbemann’s lattice $`Q_{32}`$ , \[5, Chap. 8\] has $`\tau =`$261,120.

In recent years the kissing numbers of a few other lattices in dimensions $`>32`$ have been determined. Nebe  shows that the Mordell-Weil lattice $`MW_{44}`$ has $`\tau =`$2,708,112. Nebe  shows that a 64-dimensional lattice constructed in  is extremal 3-modular, and so by modular form theory has $`\tau =`$138,458,880. Bachoc and Nebe  give an 80-dimensional lattice with $`\tau =`$1,250,172,000. Elkies  calculated the kissing number of his lattice $`MW_{128}`$: it is 218,044,170,240, over 170 times that of $`BW_{128}`$.

In the present note, we show that an elementary construction using binary codes gives better values than all of these. However, our packings are just local arrangements of spheres around the origin: we do not know if they can be modified to produce dense infinite packings.

## 2 The construction

Let $`𝒞(n,d)`$ (resp. $`𝒞(n,d,w)`$) denote a set of binary vectors of length $`n`$ and Hamming distance $`d`$ apart (resp. and with constant weight $`w`$). The maximal size of such a set is denoted by $`A(n,d)`$ (resp. $`A(n,d,w)`$) , .

One way to achieve the kissing number $`\tau _8=240`$ in eight dimensions is to take as centers of spheres the vectors of shape $`\pm 1^8`$, with a unique support (a code $`𝒞(8,8,8)!`$) and signs taken from a $`𝒞(8,2)`$, together with the vectors of shape $`\pm 2^20^6`$, where the supports are taken from a $`𝒞(8,2,2)`$ and the signs from a $`𝒞(2,1)`$. Taking all these codes to be as large as possible, we obtain a total of

$$A(8,8,8)A(8,2)+A(8,2,2)A(2,1)=12^7+\left(\genfrac{}{}{0pt}{}{8}{2}\right)2^2=240$$

spheres touching the sphere at the origin.

Our construction generalizes this as follows. For a given dimension $`n`$, we choose a sequence of support sizes $`n_0`$, $`n_1,\mathrm{},n_\mu `$ satisfying

$$nn_04n_14^2n_2\mathrm{}4^\mu n_\mu 1.$$
(1)
The $`\nu ^{\mathrm{th}}`$ set of centers that we use, for $`0\nu \mu `$, consists of vectors of shape $`\pm a_\nu ^{n_\nu }0^{nn_\nu }`$, where $`a_\nu =\sqrt{n_0/n_\nu }`$, the supports are taken from a $`𝒞(n,n_\nu ,n_\nu )`$ and the signs from a $`𝒞(n_\nu ,\frac{n_\nu }{4})`$. With optimal codes, the total number of vectors is

$$\underset{\nu =0}{\overset{\mu }{}}A(n,n_\nu ,n_\nu )A(n_\nu ,\frac{n_\nu }{4}).$$
(2)
It is easy to check that all vectors have length $`\sqrt{n_0}`$, and that by (1) the distance between any two distinct vectors is $`\sqrt{n_0}`$. It follows that (2) is a lower bound on $`\tau _n`$.

### Remarks

(1) Even if we do not know the exact values of $`A(n,d,w)`$ and $`A(n,d)`$ mentioned in (2), we can replace them by any available lower bounds, and still obtain a lower bound on the kissing number $`\tau _n`$. There is some freedom in choosing the $`n_\nu `$, which helps to compensate for our ignorance.

(2) A table of lower bounds on $`A(n,d)`$ has been given by Litsyn , extending the table in . A table of lower bounds on $`A(n,d,w)`$ for $`n28`$ is given in , but for larger $`n`$ little is known. A very incomplete table for $`n>28`$ can be found in .

(3) The construction gives a set of points on a sphere with angular separation of $`60^{}`$. It can obviously be modified to produce spherical codes with other angles.

## 3 Examples

We illustrate the construction by giving new records in dimensions 32, 36, 40, 44, 64, 80 and 128. For other examples see , and for further details about the codes see , .

#### n=32.

We take $`n_0=32`$, $`n_1=8`$, $`n_2=2`$ and use $`A(32,8)2^{17}`$ from , $`A(32,8,8)1117`$ from the complement of a lexicographic code $`𝒞(32,8,24)`$ (cf. ), obtaining a kissing number of $`A(32,32,32)A(32,8)+A(32,8,8)A(8,2)+A(32,2,2)A(2,1)12^{17}+11172^7+\left(\genfrac{}{}{0pt}{}{32}{2}\right)2^2=`$276,032.

#### n=36.

Let the 36 coordinates be labeled $`(i,j)`$, $`0i,j5`$, and let the symmetric group $`S_6`$ act by $`(i,j)(i^\pi ,j^{\sigma (\pi )})`$, where $`\pi S_6`$ and $`\sigma `$ is the outer automorphism of $`S_6`$. One can find a set of 17 orbits under the alternating group $`A_6`$, of sizes ranging from 45 to 360, whose union forms a constant weight code showing that $`A(36,8,8)2385`$. We take $`n_0=32`$, $`n_1=8`$, $`n_2=2`$ and obtain a kissing number of $`A(36,32,32)A(32,8)+A(36,8,8)A(8,2)+A(36,2,2)A(2,1)12^{17}+23852^7+\left(\genfrac{}{}{0pt}{}{36}{2}\right)4=`$438,872.

An alternative approach can be based on Warren D. Smith’s discovery (personal communication, May 1997) that the 2754 minimal vectors of the self-dual length 18 distance 8 code over $`F_4`$  yields $`\tau _{36}27542^7=`$352,512 by changing any even number of signs. By adjoining additional vectors with fractional coordinates R. H. Hardin and N. J. A. Sloane increased this to 386,570, which held the record until it was overtaken by the present construction. It is quite possible that with better clique-finding the $`F_4`$ approach will regain the lead.

#### n=40.

We take $`n_0=40`$, $`n_1=8`$, $`n_2=2`$, use a lexicographic code for $`A(40,8,8)`$, and obtain $`A(40,40,40)A(40,10)+A(40,8,8)A(8,2)+A(40,2,2)A(2,1)1589824+31162^7+\left(\genfrac{}{}{0pt}{}{40}{2}\right)2^2=`$991,792.

#### n=44.

$`A(44,44,44)A(44,11)+A(44,8,8)A(8,2)+A(44,2,2)A(2,1)12^{21}+66222^7+\left(\genfrac{}{}{0pt}{}{44}{2}\right)4=`$2,948,552.

#### n=48.

In 48 dimensions the three known extremal unimodular lattices ,  have kissing number 52,416,000. Our present construction gives less than half this value.

#### n=64.

The words of weight 16 in an extended cyclic code $`𝒞(64,16)`$ of size $`2^{28}`$ from  show that $`A(64,16,16)`$30,828. In this way we obtain a kissing number of 331,737,984.

#### n=80.

By taking 4 orbits under $`L_2(79)`$ we obtain $`A(80,16,16)`$143,780. We take $`n_0=64`$, $`n_1=16`$, $`n_2=4`$, $`n_3=1`$ and obtain $`\tau `$1,368,532,064.

#### n=128.

This is the most dramatic improvement, so we give a little more detail. Our construction uses:

$$\begin{array}{ccc}A(128,128,128)A(128,32)\mathrm{vectors}\pm 1^{128}:\hfill & \hfill & 12^{43}\hfill \\ A(128,32,32)A(32,8)\mathrm{vectors}\pm 2^{32}0^{96}:\hfill & \hfill & 5120642^{17}\hfill \\ A(128,8,8)A(8,2)\mathrm{vectors}\pm 4^80^{120}:\hfill & \hfill & 27045922^7\hfill \\ A(128,2,2)A(2,1)\mathrm{vectors}\pm 8^20^{126}:\hfill & \hfill & \left(\genfrac{}{}{0pt}{}{128}{2}\right)4\hfill \\ & & \\ \multicolumn{2}{c}{\text{for a total of}}& \text{8,863,556,495,104}\hfill \end{array}$$

Here $`A(128,32)2^{43}`$ comes from a BCH code \[10, p. 267\], $`A(128,32,32)`$512064 from a union of two orbits under $`L_2(127)`$, $`A(32,8)2^{17}`$ from , and $`A(128,8,8)`$2704592 is obtained by shortening a $`𝒞(129,8,8)`$ of size 2883408 formed from the union of 11 orbits of size 262128 under $`L_2(128)`$. The result is more than 40 times that of the Mordell-Weil lattice.

We do not expect any of these new records to survive for long, since our lower bounds for $`A(n,d)`$ and $`A(n,d,w)`$ are very weak. However, it is interesting that such a simple construction gives such dramatic improvements over the kissing numbers of the best lattices known.

Postscript. Victor Zinoviev has pointed out to us that in 1992 he and T. Ericson  obtained kissing numbers of 858800 in dimension 40 and 273935235 in dimension 64 by similar methods. Warren D. Smith (unpublished) used these methods to construct spherical codes in the 1980’s.


format: text

text:
# 1 Introduction

## 1 Introduction

Photosynthesis is the main source of energy in our biosphere. Although light-harvesting systems show a great variety in their design and function, most of them are membrane proteins comprised of a network of pigment antennae surrounding a reaction center . The antenna complex is responsible for the absorption of light and the transfer of the resulting electronic excitation to a so-called special pair of chlorophylls in the reaction center, where a charge separation across the cell membrane is initiated. This charge separation is later utilized by the cell to store energy through synthesis of ATP in a more stable form.

Oxygenic photosynthetic species employ two transmembrane protein-pigment complexes for light-harvesting, named photosystems I and II. Photosystem I (PSI) is a ubiquitous protein-pigment complex found in green plants, algae and cyanobacteria. Located in the bacterial membrane, it absorbs sun light and uses its energy to transfer electrons across the cell membrane. PSI is remarkable in that it contains the antenna complex, reaction center and the electron transfer chain all within the same protein. For a recent review on PSI physiology and function we refer the reader to Refs. .

A high resolution atomic structure was recently obtained for PSI from the cyanobacterium Synechococcus (S.) elongatus . The structure exhibits an aggregate of 96 chlorophyll molecules, including the special pair, P700, at which the primary charge separation takes place. This finding opens up the opportunity to study in detail the light-harvesting function of PSI, i.e., one can computationally examine the light-harvesting pathway from initial absorption of a photon to the transfer of electrons by the special pair. In the past, a similar leap in structural information on the peripheral light-harvesting complexes LH1 and LH2 of purple bacteria resulted in the photosynthetic unit of purple bacteria being the best understood antenna complex today .

It is of interest to compare the structure and function of PSI to the photosynthetic unit of purple bacteria. Both systems contain on the order of a hundred chromophores per reaction center. In purple bacteria these chromophores are organized in highly symmetric ring-like structures . In contrast, the chlorophyll aggregate in PSI forms a rather random array surrounding the reaction center (see Fig. 1), the only hint of order being a pseudo-$`C_2`$ symmetry apparent in the structure.

Purple bacteria have the oldest known photosynthetic apparatus, while cyanobacteria and their relatives are evolutionarily more recent . It is of interest to compare the light-harvesting systems of these lifeforms and to elucidate the drive toward a more efficient system in the more recently developed light-harvesting complexes.

A related issue is whether the peculiar arrangement of chlorophylls found in the reported structure is essential for the functioning of PSI or whether the system can be described successfully assuming some degree of randomness. For example, is every single chlorophyll essential for the proper functioning of the system? This question also constitutes a modeling challenge. Due to its circular symmetry and repetitive elements a purple bacterial light-harvesting complex can be modeled with only a few physical parameters describing chlorophyll site energies and electronic couplings . PSI, on the other hand, features 96 chlorophylls that show few regularities and seemingly require a large number of parameters for their description.

There are two possible ways to approach the complexity of the chlorophyll aggregate in PSI. In one approach, typically taken in chemistry, one strives for an ever increasing accuracy in a model through a large number of specific parameters. This approach implies that the specific properties of every one of the 96 chlorophylls in PSI matter. In another approach, typically taken in physics, one identifies a minimal set of characteristics, e.g., average properties of the 96 chlorophylls, relevant for the description of function and observed spectra.

The most relevant functional property of PSI is its high efficiency in using excitation energy for charge separation in the reaction center despite thermal disorder and ever present internal damage. One seeks to know how robust PSI is against perturbations such as ‘pruning’ or ‘bleaching’ of individual chlorophylls or changes in the protein that effect chlorophyll site energies. Robustness in biology is not a new concept. The issue was addressed earlier in Refs.  and, recently in the study of biochemical networks involved in bacterial chemotaxis . These studies suggest that biological networks exhibit a high degree of insensitivity to fine-tuning of physical parameters.

Another issue of interest is the optimality of the design of the chlorophyll aggregate in PSI. Is it possible to rearrange or reorient the pigments in PSI to achieve a better quantum yield or did evolution establish an optimal arrangement? This is likely easier to answer in silico than in vivo.

While the issues of robustness and optimality can be addressed well by the physicists’ approach, details of the low temperature excitation transfer dynamics and spectral features are likely better studied by the chemists’ approach. The longest wavelength features of the spectrum attributed to so-called ‘red chlorophylls’ are an example of this. Red chlorophylls are responsible for broadening the absorption profile of PSI, but their role is not fully understood. For a review we refer the reader to Ref. .

In order to answer the questions raised above we study ensembles of various chlorophyll arrangements. This permits us to investigate the dependence of the light harvesting function of PSI in detailed or coarse descriptions as well as the robustness and optimality of the system. For this purpose, we use the high resolution structure reported in  as a starting point to describe the spectral and excitation transfer characteristics of PSI. The structure reveals the relative orientation of chlorophyll molecules, which makes it possible to construct an effective Hamiltonian for the chlorophyll aggregate. The resulting quantum yield for charge separation and the average excitation lifetime is determined based on this Hamiltonian.

We also investigate the excitation transfer dynamics that results after the absorption of a photon by PSI. Besides determining the overall lifetime of an excitation and the quantum yield (the probability that photon absorption leads to electron transfer from the special pair), one wishes to know how excitation migrating through PSI repeatedly reaches and escapes the special pair, i.e., engages in sojourns from the special pair into the chlorophyll pool. We suggest a kinetic expansion method, the sojourn expansion, to account for the outlined excitation migration in a light-harvesting complex.

The organization of this paper is as follows: In the next section we introduce the effective Hamiltonian for the chlorophyll aggregate of PSI. Sec. 3 utilizes this Hamiltonian to discuss candidates for the red chlorophyll states. In Sec. 4 we discuss the effect of thermal disorder on the spectral properties of PSI in terms of random matrix theory. In Sec. 5 we introduce the general theoretical framework for the study of excitation transfer dynamics, including the sojourn expansion, and apply it to PSI. The robustness of the efficiency of the system is investigated under the effects of the fluctuations in chlorophyll site energies and the pruning of individual chlorophylls. In Sec. 6 we examine the optimality of the arrangement of chlorophylls in PSI. Sec. 7 contains a general discussion and our conclusions. In an appendix we summarize our method to compute inter-chlorophyll couplings.

## 2 Effective Hamiltonian for chlorophyll aggregate

In this section we construct an effective Hamiltonian for the aggregate of chlorophyll molecules in PSI defined in a basis of single chlorophyll excitations. Isolated cyanobacterial PSI is known to exist in vivo mostly as a trimer . In this paper we will examine the behavior of the monomeric complex. Of the 96 chlorophylls present in PSI, 95 are chlorophyll $`a`$ molecules and one of the special pair chlorophylls is a chlorophyll $`a^{}`$. For computational simplicity, we will assume below all chlorophylls to be of chlorophyll $`a`$ type.

The lowest excited state of a chlorophyll molecule is the $`Q_y`$ state . A set of basis states for an effective Hamiltonian is then given by

$$|i=|\varphi _1\varphi _2\mathrm{}\varphi _i^{}\mathrm{}\varphi _N,i=1,2,\mathrm{}N,$$
(1)
in which the $`i`$th chlorophyll is $`Q_y`$-excited, all the other chlorophylls being in their ground states. The total number of chlorophylls is $`N=96`$. By using this basis, an effective Hamiltonian can be written as

$`H_0=\left(\begin{array}{ccccc}ϵ_1& W_{12}& W_{13}& \mathrm{}& W_{1N}\\ W_{21}& ϵ_2& W_{23}& \mathrm{}& W_{2N}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ W_{N1}& W_{N\mathrm{\hspace{0.17em}2}}& W_{N\mathrm{\hspace{0.17em}3}}& \mathrm{}& ϵ_N\end{array}\right).`$ (6)
Here, $`ϵ_i`$ denotes the chlorophyll site energy for the $`Q_y`$ state of the $`i`$th chlorophyll, which needs to be determined later, and $`W_{ij}`$ denotes the coupling between the $`i`$th and the $`j`$th chlorophylls.

The electronic coupling, $`W_{ij}`$, between excited states of two chlorophyll molecules has two contributions

$$W_{ij}=W_{ij}^c+W_{ij}^{ex},$$
(7)
where $`W_{ij}^c`$ corresponds to a direct Coulomb term  and $`W_{ij}^{ex}`$ corresponds to an electron exchange term . In PSI, almost all of the chlorophylls are at a Mg-Mg distance of greater than 7 Å, a range in which Coulomb interaction dominates ; therefore we keep only the contribution of the direct Coulomb term, $`W_{ij}^c`$, in the remainder of the paper.

The electronic couplings due to the direct Coulomb interaction are calculated through two approaches. In the first (approximate) approach, the couplings are determined assuming a dipole-dipole interaction between transition dipole moments of the chlorophyll $`Q_y`$ states. In this approximation, the couplings are given by (see, for example, Ref. )

$`W_{ij}=C\left({\displaystyle \frac{𝐝_i𝐝_j}{r_{ij}^3}}{\displaystyle \frac{3\left(𝐫_{ij}𝐝_i\right)\left(𝐫_{ij}𝐝_j\right)}{r_{ij}^5}}\right)`$ (8)
where $`𝐝_i`$ are unit vectors along transition dipole moments from the ground state to the $`Q_y`$ state of the $`i`$th chlorophyll and $`𝐫_{ij}`$ is the vector connecting the Mg atoms of chlorophylls $`i`$ and $`j`$; the transition dipole moments for the $`Q_y`$ state are taken to be along the vector connecting the NB and ND atoms. The latter assumption is seen to be valid by comparing the transition dipole moment unit vectors thus computed to the ones obtained in the course of the full Coulomb semi-empirical method introduced below. The atomic positions are taken from crystallographic data . The parameter $`C=116000\mathrm{\AA }^3\mathrm{cm}^1`$ in Eq. (8) is determined by enforcing asymptotic agreement at large distances with the full Coulomb couplings introduced below . Full Coulomb couplings in turn are fixed by the oscillator strength of the $`Q_y`$ transition, which is taken to be $`23`$ Debye2 . Since the transition dipole moment vectors are computed simply from geometry, application of Eq. (8) has a low computational cost. Naturally, the dipolar approximation is the more accurate, the larger the distance between the chlorophylls.

The second approach, following Ref. , computes electronic couplings using a semi-empirical Pariser-Parr-Pople (PPP) Hamiltonian as utilized in Refs. . In this approach, the contribution of all orders of multipoles are included, i.e., not only the dipole-dipole contribution. Therefore, we will refer to this description as the full Coulomb (FC) approach. The details of this approach are outlined in the appendix.

In Fig. 2 are illustrated the first five eigenstates of the effective Hamiltonian (6) with full Coulomb couplings and identical site energies. Excitons are localized typically only over a few chlorophylls, in some cases coinciding with strongly coupled pairs.

## 3 Red chlorophyll candidates for photosystem I

A number of chlorophyll molecules in PSI of S. elongatus are known to absorb light at longer wavelengths than the special pair P700. The number of these ‘red chlorophylls’ is suggested to be between seven and eleven . In Ref.  individual absorption peaks for the red chlorophyll states were observed for 4 K at wavelengths 708 nm, 715 nm, and 719 nm. Among possible causes for this red shift in the absorption peak of these chlorophylls are the effect of the local protein environment on the site energies and the effect of excitonic splitting due to a strong coupling between pairs of chlorophylls.

In Ref.  candidates for the red chlorophylls were suggested based on the orientation of chlorophylls and the assumption of a strong excitonic splitting as the cause of the spectral red shift. In Ref.  a match between these candidates and the experimentally resolved red absorption peaks were suggested. In this section we extend this discussion by comparing the dipolar couplings with the corresponding full Coulomb couplings. In Ref.  a trimer (B31-B32-B33) and three dimers (A32-B7, B37-B38, A38-A39) of chlorophylls were given as red chlorophyll candidates (see Fig. 3; nomenclature of Ref.  is used for naming chlorophylls). Indeed in the framework of the dipolar couplings these chlorophylls form the strongest coupled pairs. However, in the framework of the full Coulomb couplings a different picture emerges. This may not be surprising as the two approaches tend to disagree over short distances, such as the inter-chlorophyll distances within strongly coupled pairs.

In the full Coulomb description the couplings of the chlorophyll pairs in the trimer B31-B32-B33 are significantly (more than a factor 3) smaller than their corresponding dipolar couplings. Since there are chlorophyll pairs with much stronger couplings in PSI, the trimer loses its position as a red chlorophyll candidate in the full Coulomb picture.

Furthermore, a strong coupling is not enough to cause a red shift in the absorption spectrum. Viewed as a two state system, a chlorophyll pair will have one higher energy excitonic state and one lower energy excitonic state, each with a certain oscillator strength. A red shift can be associated only with a case where the lower energy excitonic state has a higher oscillator strength than the higher energy excitonic state. Otherwise, there needs to be a significant shift in the chlorophyll site energies due to the protein environment to account for the red shift. It turns out that the two strongest coupled pairs, B37-B38 and A38-A39, have most of their oscillator strength in the higher energy state. Therefore, we leave them out of the list of red chlorophyll candidates.

Four pairs of chlorophylls, which have the strongest couplings as well as the higher oscillator strength in the lower lying state, are suggested as red chlorophyll candidates (highlighted in Table 1 and Fig. 3). A match between the observed red absorption peaks and the suggested candidates cannot be done reliably without knowledge of the site energies. However, matching the furthest red shift with the strongest coupling we tentatively suggest that A32-B7 is responsible for the 719 nm absorption peak, A33-A34 and A24-A35 for the 715 nm peak, and B22-B34 for the 708 nm peak.

It must be noted that excitonic splitting alone cannot be the only cause for observed red shifts. For example, in case of two identical chlorophylls, which have their $`Q_y`$ absorption peaks at 675 nm, a coupling of about 900 cm-1 would be needed to place the lower energy excitonic state at 719 nm. This suggests that the red chlorophyll states are caused by the contribution of both excitonic splitting and site energy shifts due to the local protein environment. Furthermore, any isolated chlorophyll may contribute to the red chlorophyll band of the PSI spectrum solely due to its site energy shift, a possibility which is totally neglected in our present picture.

The trimer of chlorophylls, B31-B32-B33, deserves a closer inspection than presented thus far. As seen in Fig. 4, the atomic structure indicates that the three chlorophylls are connected to each other by water molecules, suggesting that they might have to be treated computationally as a supermolecule. In regard to the possible failure of our description, we note that the His residue coordinating chlorophylls B31-B32-B33 is absent in Synechocystis, which exhibits a weaker long wavelength absorption band. Indeed, the close packing of chlorophylls B31-B32-B33 is reminiscent of J-aggregates , which are known to have red shifted spectra compared to the spectra of their building blocks. We may suggest, hence, that the chlorophylls B31-B32-B33 may yet be responsible for red absorption, despite the weak couplings emerging from our description.

## 4 Effect of thermal disorder on spectral properties

Spectral properties and excitation transfer dynamics of a photosynthetic system are two complementary sources of information. Both are determined by the Hamiltonian describing the system. However, the spectral properties are influenced more by chlorophyll site energies, i.e., the diagonal entries of the effective Hamiltonian, whereas the excitation transfer dynamics more directly reflects the role of inter-chlorophyll couplings.

A reconstruction of the (low temperature) absorption spectrum from the Hamiltonian is not possible without a knowledge of the site energies. This is seen, for example, in the much smaller seven chlorophyll system of the FMO complex . Similarly the inverse problem of using the spectrum to reconstruct the site energies (assuming the couplings are known) is ill posed, since with as many parameters as the number of chlorophylls in the system one cannot unambiguously assign the site energies. Indeed it is a very important test for any calculation of the chlorophyll site energies to faithfully reproduce the experimentally observed absorption spectrum, especially at low temperatures.

Nevertheless spectral data can provide valuable information on the general distribution of chlorophyll site energies if not on individual sites. Therefore, we study the spectral properties of PSI in the light of the effective Hamiltonian introduced in Sec. 2, before we proceed with a detailed discussion of excitation transfer dynamics.

### 4.1 Random matrix theory as a description of thermal disorder

Both the spectral features and the excitation transfer dynamics in a photosynthetic system need to be described in the presence of strong thermal disorder. Standard zero-temperature quantum theory remains inadequate in this respect. Below we present a method, based on random matrix theory , to study static disorder in a system described by an effective Hamiltonian, such as the one given in Sec. 2. For a more detailed account of the application of random matrix theory to the description of static disorder in photosynthetic systems the reader is referred to Ref. . For a discussion of dynamic disorder in a light-harvesting complex the reader is referred to Refs. .

We consider an ensemble of Hamiltonians that describes a set of photosynthetic systems in thermal equilibrium. In contrast to a canonical ensemble from statistical mechanics, whose members are points in the phase space of one Hamiltonian, we consider an ensemble of Hamiltonians each of which represents another copy of the photosynthetic system in question. We will be interested in spectral properties as averaged over this ensemble.

For a description of this ensemble, we employ the sum of the noise-free effective Hamiltonian, $`H_0`$, based on the crystal structure, and a random part, $`R`$, representing thermal fluctuations

$`H=H_0+R.`$ (9)
Here the matrix $`R`$ is drawn from a certain probability distribution, $`P\left(R\right)`$, representing the effects of thermal disorder. All relevant spectral quantities are defined in terms of an ensemble average over the random part, $`R`$, with respect to the weight function, $`P\left(R\right)`$.

The density of states and the directionally averaged absorption spectrum for (9) are defined by

$`\rho \left(\omega \right)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left(\omega E_i\right),`$ (10)
$`\alpha \left(\omega \right)`$ $`=`$ $`{\displaystyle \frac{4\pi ^2\omega n}{3c}}{\displaystyle \underset{i}{}}\left|𝐃_i^2\right|\delta \left(\omega E_i\right),`$ (11)
where $`<\mathrm{}>d\left[R\right]P\left(R\right)\mathrm{}`$ and $`E_i`$ are the eigenvalues of $`H=H_0+R`$. The transition dipole moments for eigenstates, $`𝐃_i=_mc_i\left(m\right)𝐝_m`$, are given in terms of the transition dipole moment unit vectors of individual sites, $`𝐝_m`$, and the expansion coefficients, $`c_i\left(m\right)`$, of eigenstates.

These definitions for the average spectral behavior are superior to the ‘Gaussians on sticks’ approach sometimes utilized, where Gaussian envelopes are put around the eigenvalues of $`H_0`$ to reconstruct the spectrum. Especially in cases where there is a significant overlap between the envelopes of multiple eigenvalues, the average spectral properties are seen to deviate significantly. As an extreme example one can consider a case of $`N`$ degenerate eigenvalues . In this case a sum of Gaussian envelopes gives another Gaussian with the same width, whereas the width of the average spectrum should be proportional to $`\sqrt{N}`$.

The nature of the probability distribution, $`P\left(R\right)`$, as a function of temperature is not easy to determine. However, studies in random matrix theory have shown that many average spectral features are largely independent of the exact shape of $`P\left(R\right)`$ and instead depend mainly on its width . This kind of independence from the specific probability distribution is akin to the central limit theorem and is known as random matrix universality. This universality allows one to make simplifying assumptions on the nature of $`P\left(R\right)`$ .

Although most results on random matrix universality are based on the limit where the matrix size becomes large, remarkable spectral similarity persists even for mesoscopically sized ensembles. In an earlier investigation, we have applied random matrix universality to the peripheral light-harvesting complex, LH2, of purple bacteria , which has 16 strongly coupled bacteriochlorophylls. It was seen that, when the width of the disorder term, $`P\left(R\right)`$, is taken into account the average density of states and the absorption spectrum are rather insensitive to the changes of $`P\left(R\right)`$. Analytical formulae for the density of states at finite temperature as a function of the zero temperature spectrum were derived for a simple form of $`P\left(R\right)`$ in Ref. .

In the following we take the probability distribution, $`P\left(R\right)`$, to describe diagonal disorder

$`P_{diag}\left(R\right)`$ $`=`$ $`N_{diag}{\displaystyle \underset{i}{}}\mathrm{exp}\left({\displaystyle \frac{1}{2\nu _{diag}^2}}R_{ii}^2\right){\displaystyle \underset{ij}{}}\delta \left(R_{ij}\right).`$ (12)
For the width of the distribution (12) we take $`\nu _{diag}=70`$ cm-1 at 4 K . A similar description has been used in Ref.  for the light harvesting system of purple bacteria. We use the parameter $`\nu _{diag}`$ to numerically evaluate the absorption spectrum as defined in (11). In Fig. 5 is shown the absorption spectrum at 4 K based on the effective Hamiltonian with full Coulomb couplings and homogeneous site energies corresponding to 675 nm. Ignoring the red chlorophyll band of the spectrum, a width for the bulk of about 230 cm-1 (corresponding to a FWHM of 550 cm-1) at 4 K is reported in Ref. . However, the width for the 4 K spectrum in Fig. 5 measures only about 140 cm-1. The mismatch can be attributed to our neglect of the site energy heterogeneity.

Though it is not possible to reconstruct the chlorophyll site energies from the spectrum, one can obtain an estimate on the width of the distribution of site energies. In order to accomplish this we will ignore for the moment the red chlorophyll band of the spectrum and instead concentrate on the main peak as given in Ref. . The width of the absorption spectrum for the actual (inhomogeneous) site energies, $`W_{inhom}`$, the width of the absorption spectrum for the homogeneous site energies, $`W_{hom}`$, and the width of the distribution of site energies (heterogeneity), $`W_ϵ`$, are approximately related to each other by

$`W_{inhom}^2=W_{hom}^2+W_ϵ^2.`$ (13)
As discussed above, assuming $`W_{inhom}=230`$ cm-1 and $`W_{hom}=140`$ cm-1 results in a site energy distribution width of about $`W_ϵ=180`$ cm-1 (FWHM of 430 cm-1). The validity of this assertion can be numerically verified by generating an ensemble of Hamiltonians with the same full Coulomb couplings but with random site energies for the given value of $`W_ϵ`$. This value will be utilized later in the discussion of excitation transfer dynamics.

## 5 Excitation transfer dynamics in photosystem I

In this section we construct the excitation transfer rates between chlorophylls using Förster theory . We will follow the methodology suggested in Ref. . The excitation transfer rate from chlorophyll $`i`$ to chlorophyll $`j`$ is

$`T_{ij}={\displaystyle \frac{2\pi }{\mathrm{}}}\left|H_{ij}\right|^2J_{ij},J_{ij}={\displaystyle S_i^D\left(E\right)S_j^A\left(E\right)𝑑E},`$ (14)
where $`H_{ij}`$ is the coupling between the chlorophylls and $`J_{ij}`$ is the spectral overlap between the emission spectrum, $`S_i^D\left(E\right)`$, of the donor chlorophyll $`i`$ and the absorption spectrum, $`S_j^A\left(E\right)`$, of the acceptor chlorophyll $`j`$. Following Refs.  we approximate these spectra by Gaussians

$`S_i^D\left(E\right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\nu }}\mathrm{exp}\left({\displaystyle \frac{\left(E_iSE\right)^2}{2\nu ^2}}\right),`$
$`S_j^A\left(E\right)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }\nu }}\mathrm{exp}\left({\displaystyle \frac{\left(E_jE\right)^2}{2\nu ^2}}\right),`$ (15)
where $`E_i`$ and $`E_j`$ are the absorption peaks for the chlorophylls and $`S`$ is the Stokes shift. We take the Stokes shift to be equal to 160 cm-1 at room temperature  and 20 cm-1 at 4 K . We have assumed identical widths for the emission and absorption spectra of 240 cm-1 at room temperature . In the case of identical site energies, rates for forward transfer and back transfer are equal, $`T_{ij}=T_{ji}`$, but this is not the case when the site energies differ. The network of excitation transfer rates between chlorophylls of PSI for identical site energies are shown in Fig. 6 for both the full Coulomb and dipolar Hamiltonians.

The largest transfer rate away from a given chlorophyll suggests the average excitation lifetime at that chlorophyll. Fig. 7 shows the largest excitation transfer rates for all 96 chlorophylls. It can be seen that these rates correspond to an average largest transfer rate between chlorophylls of about 11 ps-1.

Using the transfer rates between the chlorophylls, a master equation can be constructed for the rate of change of probabilities, $`p_i\left(t\right)`$, describing the likelihood that the chlorophyll $`i`$ is electronically excited at time t. The corresponding rate equation is

$`{\displaystyle \frac{d}{dt}}p_i\left(t\right)`$ $`=`$ $`{\displaystyle \underset{j}{}}K_{ij}p_j\left(t\right),`$ (16)
$`K_{ij}`$ $`=`$ $`T_{ji}\delta _{ij}\left({\displaystyle \underset{k}{}}T_{ik}+k_{diss}+\delta _{i,P700}k_{CS}\right),`$ (17)
where $`k_{diss}`$ denotes the dissipation (internal conversion) rate, which is assumed to be uniform among the chlorophylls, $`k_{CS}`$ denotes the charge separation rate at the special pair P700 and $`\delta _{i,P700}`$ is equal to 1 when $`i`$ is one of the two chlorophylls in the special pair and zero otherwise. The dissipation and charge separation rates are not known to great accuracy. Inspired by the purple bacterial light-harvesting systems  we assume a dissipation rate of $`k_{diss}=\left(1\mathrm{ns}\right)^1`$. The charge separation rate at P700 is observed to be between 1 ps-1 and (3 ps)-1 . In the following we will assume a charge separation rate of $`k_{CS}=\left(1.5\mathrm{ps}\right)^1`$.

Given an initial distribution, $`p_i\left(0\right)`$, the formal solution to Eq. (16) can be written

$`|p\left(t\right)=\mathrm{exp}\left(Kt\right)|p\left(0\right).`$ (18)
Using this expression, explicit formulae for the average excitation lifetime, the quantum yield, and overall dissipation rate can be given .

Let us define a uniform (non-normalized) state, $`|\mathrm{𝟏}>_i|i>`$, describing equally likely occupation probabilities for all chlorophylls. If we assume this uniform distribution for the initial state, i.e., $`|p\left(0\right)=N^1|\mathrm{𝟏}>`$, then the average excitation lifetime, $`\tau `$, the quantum yield, $`q`$, and the overall dissipation probability, $`d`$, are all given in terms of similar expressions

$`\tau `$ $`=`$ $`{\displaystyle \frac{1}{N}}\mathrm{𝟏}\left|K^1\right|\mathrm{𝟏},`$ (19)
$`q`$ $`=`$ $`{\displaystyle \frac{1}{N}}k_{CS}P700\left|K^1\right|\mathrm{𝟏}=1d,`$ (20)
$`d`$ $`=`$ $`{\displaystyle \frac{1}{N}}k_{diss}\mathrm{𝟏}\left|K^1\right|\mathrm{𝟏}=k_{diss}\tau ,`$ (21)
where $`|P700_i\delta _{i,P700}|i`$.

Average excitation lifetimes and the quantum yields for various models are provided in Table 2. The data in the first three rows are based on an effective Hamiltonian that assumes identical site energies for the two chlorophylls of the P700 pair. The site energies were chosen such that the P700 excitonic state with the highest oscillator strength coincides with the 698 nm (14327 cm-1) absorption peak at 4 K. Since one of the two P700 chlorophylls is a chlorophyll $`a^{}`$, the identical site energy assumption for the P700 is likely to be unrealistic, however it shall be adequate for an approximate description of the excitation transfer dynamics. For full Coulomb couplings the highest oscillator strength lies at the higher energy excitonic state (see Table 1); since the respective coupling between the two P700 chlorophylls is 47.6 cm-1, the P700 site energy is chosen to be 14279 cm-1. For dipolar couplings the highest oscillator strength lies at the lower energy excitonic state (the sign of the coupling is different for the two methods); the coupling between the two P700 chlorophylls is 272 cm-1 requiring a P700 site energy of 14599 cm-1. The site energies for the remaining chlorophylls are placed at 675 nm (14815 cm-1) unless otherwise noted.

The efficiency of 97.3 % and the average excitation lifetime of 27.4 ps for full Coulomb couplings, given in the first row of Table 2, are in general agreement with the near unit efficiency and the 20-40 ps lifetime reported for PSI at room temperature . A recent estimate of the average excitation lifetime reported for trimeric PSI in S. elongatus is 35.8 ps . A comparison between the first and the last rows in Table 1 shows the importance of having P700 chlorophylls at a lower energy than the bulk of the chlorophylls in order to decrease the trapping time and to increase efficiency.

In Fig. 8 are provided the distribution of quantum yield and average excitation lifetime across an ensemble of Hamiltonians with full Coulomb couplings, P700 tuned to the observed absorption peak, and all other chlorophyll site energies randomly distributed with a width of 180 cm-1 around 675 nm (thus the ensemble covers different realizations of heterogeneity, not static disorder). Fig. 8 reveals that, at room temperature, fluctuations in the site energies (heterogeneity) have no considerable effect on efficiency. This allows us to reconstruct room temperature excitation transfer dynamics with reasonable accuracy even without detailed knowledge of the chlorophyll site energies. However, at 4 K, the line shapes of individual chlorophylls become so narrow that the spectral overlap integrals in Eq. (14) become unrealistically small for misplaced site energies. Also at low temperatures the broadening of spectral line shapes due to electron-phonon couplings  need to be taken into account as this broadening is likely to effect the overlap integrals in the context of Förster theory. We find it impossible to construct a reliable picture of the low temperature excitation transfer dynamics with the current data.

It is of interest to know in how far low energy chlorophylls control the excitation transfer dynamics and trapping in PSI. In this regard, one may construct models in which the red chlorophyll candidates given in Sec. 3 are used to refine the assignment of chlorophyll site energies. However, it must be kept in mind that an accurate assignment of an individual chlorophyll pair to an observed red chlorophyll state cannot be done reliably without a precise knowledge of site energies. It is seen that the room temperature excitation transfer dynamics is not effected greatly by the inclusion of this extra information in a model where all other site energies are unknown or taken to be random. Fig. 8 supports this claim. As a result we find it rather unfruitful to proliferate the multitude of models already given in Table 2. At room temperature the role of low energy chlorophylls appear to be simply to extend the absorption profile to longer wavelengths rather than to have a profound effect on the excitation transfer dynamics. However, at low temperatures red chlorophylls may function as effective traps and could serve as probes to investigate the excitation transfer dynamics of PSI in detail.

### 5.1 Effect of removal of individual chlorophylls on excitation transfer dynamics

One way to probe the function of individual chlorophylls in PSI as well as the degree of robustness in its design is to examine the effect of the removal of chlorophylls on the efficiency of the system. To this end we remove various chlorophylls one at a time and examine the quantum yield of the remaining system.

Fig. 9 shows the quantum yield as a function of the removed chlorophyll. The two P700 chlorophylls are not removed in this study. Although the effect of the removal of the other four reaction center chlorophylls on the quantum yield are shown, it is hard to test their effect directly, since they are also a part of the electron transfer chain. Therefore, we will concentrate on the remaining 90 chlorophylls.

It is interesting to note that, except for the four reaction center chlorophylls which are immediate neighbors to the P700, removal of a chlorophyll hardly has any effect on the quantum yield. In fact, for most of the chlorophylls a deletion results in an increase in the quantum yield. This is a consequence of the relative yield being defined with respect to the remaining chlorophylls, not with respect to the unperturbed system.

Thus it is seen that, outside the reaction center, no chlorophylls play the role of a gatekeeper for excitation transfer, whose removal would have an adverse effect on the yield. Instead the individual chlorophylls contribute to the overall cross-section of the antenna complex. Outside the electron transfer chain, the four chlorophylls with the highest impact on efficiency (A26, B24, B39 and A40) are the nearest neighbors of the six reaction center chlorophylls. These chlorophylls, especially A40 and B39, appear to link the antenna system to the reaction center and are therefore referred to as connecting chlorophylls. As an extreme example one may consider the removal of all four of these chlorophylls simultaneously. Even in this case the quantum yield reduces to only 96.85 $`\%`$ from the value of 97.25 $`\%`$ for the unperturbed complex with the same Hamiltonian; the average excitation lifetime increases from 27.4 ps to 31.5 ps upon removal. These results agree with some of the simulations reported in Ref. .

### 5.2 An expansion of the average lifetime in terms of return times to the reaction center

The average excitation lifetime, $`\tau `$, can be expanded in terms of processes describing the initial delivery of excitation to the special pair and possible subsequent returns following detrapping. In order to establish such an expansion, we rewrite the excitation transfer rates (16) in a form where the P700 is treated as a single excitonically coupled unit as opposed to two separate chlorophylls.

A two state Hamiltonian for P700 is given by

$`H_{P700}=\left(\begin{array}{cc}ϵ_1& U_{P700}\\ U_{P700}& ϵ_2\end{array}\right),`$ (24)
where $`ϵ_1`$ and $`ϵ_2`$ are the site energies for the two P700 chlorophylls and $`U_{P700}`$ denotes the coupling between them. Not knowing the values of the site energies it will be assumed below that $`ϵ_1=ϵ_2`$. However, the expansion below can be readily applied to nonidentical site energies.

Let us denote the eigenvalues of the two state system (24) by $`E_1`$ and $`E_2`$ and the corresponding eigenfunctions by

$`\mathrm{\Psi }_1^{P700}=c_{1,1}|1+c_{1,2}|2,`$
$`\mathrm{\Psi }_2^{P700}=c_{2,1}|1+c_{2,2}|2,`$ (25)
respectively. Then the coupling of any of the other chlorophylls, $`j=3,\mathrm{},N`$, to any eigenstate, $`m=1,2`$, of P700 is

$`\stackrel{~}{H}_{mj}=\stackrel{~}{H}_{jm}={\displaystyle \underset{\alpha =1}{\overset{2}{}}}c_{m,\alpha }H_{j\alpha },`$ (26)
where $`H_{j\alpha }`$ are the elements of the original effective Hamiltonian (6).

We will assume that thermal equilibration between eigenstates of P700 have been achieved before excitation transfer out of P700 occurs. This is likely to be a reasonable assumption at room temperature for full Coulomb couplings, where $`k_BT=209`$ cm-1 is larger than the interchlorophyll couplings. The expansion outlined below can also be formulated for a reaction center where multiple chlorophylls without a strong excitonic character are responsible for charge separation, making the thermal equilibration assumption invalid. However, in that case the final formulation of the expansion cannot be written as succinctly as Eq. (42) below.

The excitation transfer from P700 to any other chlorophyll is given by a Boltzmann weighted sum

$`\stackrel{~}{T}_{P700,j}`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{2}{}}}\stackrel{~}{T}_{mj},`$
$`\stackrel{~}{T}_{mj}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \frac{e^{E_m/k_BT}}{_{n=1}^2e^{E_n/k_BT}}}\left|\stackrel{~}{H}_{mj}\right|^2{\displaystyle S_m^D\left(E\right)S_j^A\left(E\right)𝑑E},`$ (27)
where we assume an identical line shape for the eigenstates of P700, chosen equal to the lineshape of all other chlorophylls in PSI. The excitation transfer rate from another chlorophyll to P700 is

$`\stackrel{~}{T}_{j,P700}={\displaystyle \underset{m=1}{\overset{2}{}}}T_{j,m},`$ (28)
where $`T_{j,m}`$ are given by Eq. (14), except that $`m=1,2`$ now denotes an eigenstate of P700 as given in Eq. (25).

A reduced $`\left(N1\right)\times \left(N1\right)`$ transfer rate matrix, $`\stackrel{~}{T}`$, can be constructed using these rates. The corresponding matrix, $`\stackrel{~}{K}`$, that enters in Eq. (16) governing the kinetics of PSI, is

$`\stackrel{~}{K}_{ij}`$ $`=`$ $`\stackrel{~}{T}_{ji}\delta _{ij}\left({\displaystyle \underset{k}{}}\stackrel{~}{T}_{ik}+k_{diss}+\delta _{i,P700}k_{CS}\right).`$ (29)
Average excitation lifetime, $`\tau `$, and quantum yield of the system are defined as in Eqs. (19) and (20), e.g.,

$`\tau `$ $`=`$ $`{\displaystyle \frac{1}{N1}}\mathrm{𝟏}\left|\stackrel{~}{K}^1\right|\mathrm{𝟏},`$ (30)
where $`|\mathrm{𝟏}=|P700+_{j=3}^N|j`$ denotes the (nonnormalized) uniform initial state (in the basis of this reduced Hamiltonian $`|P700`$ denotes the vector with a one in the first column and zeroes everywhere else).

We introduce now the expansion of the average excitation lifetime. The terms that arise in this expansion are introduced and explained in Fig. 10. For the purpose of this expansion we separate from $`\stackrel{~}{K}`$ the operator, $`\mathrm{\Delta }`$, that describes detrapping

$`\stackrel{~}{K}`$ $``$ $`\kappa +\mathrm{\Delta },`$
$`\mathrm{\Delta }`$ $``$ $`\left(\begin{array}{cccc}0& 0& \mathrm{}& 0\\ \stackrel{~}{T}_{P700,3}& 0& \mathrm{}& 0\\ \stackrel{~}{T}_{P700,4}& 0& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \stackrel{~}{T}_{P700,N}& 0& \mathrm{}& 0\end{array}\right).`$ (36)
The inverse, $`\stackrel{~}{K}^1`$, that enters in Eq. (30), can be expressed through Taylor expansion in $`\mathrm{\Delta }`$

$`\stackrel{~}{K}^1`$ $`=`$ $`\left(\kappa +\mathrm{\Delta }\right)^1=\left(\mathrm{𝟏}+\kappa ^1\mathrm{\Delta }\right)^1\kappa ^1`$ (37)
$`=`$ $`\kappa ^1\kappa ^1\mathrm{\Delta }\kappa ^1+\kappa ^1\mathrm{\Delta }\kappa ^1\mathrm{\Delta }\kappa ^1\mathrm{}`$
where $`\mathrm{𝟏}`$ is the $`\left(N1\right)\times \left(N1\right)`$ identity matrix. The validity of this expansion can be proven by noting that the absolute value of the eigenvalues of $`\kappa ^1\mathrm{\Delta }`$ are less than one. Substituting this expansion into Eq. (30) results in a decomposition of the lifetime into a series

$`\tau `$ $``$ $`\tau _0+\tau _1+\tau _2+\mathrm{},`$ (38)
where

$`\tau _0`$ $`=`$ $`{\displaystyle \frac{1}{N1}}\mathrm{𝟏}\left|\kappa ^1\right|\mathrm{𝟏},`$
$`\tau _1`$ $`=`$ $`+{\displaystyle \frac{1}{N1}}\mathrm{𝟏}\left|\kappa ^1\mathrm{\Delta }\kappa ^1\right|\mathrm{𝟏},`$
$`\tau _2`$ $`=`$ $`{\displaystyle \frac{1}{N1}}\mathrm{𝟏}\left|\kappa ^1\mathrm{\Delta }\kappa ^1\mathrm{\Delta }\kappa ^1\right|\mathrm{𝟏},`$ (39)
$`\mathrm{}`$
One can simplify Eq. (39) by noting that $`\mathrm{\Delta }`$ can be expressed as $`\mathrm{\Delta }=W_D|TP700|`$ where $`W_D=_{j=3}^NT_{P700,j}`$ is the total detrapping rate from the special pair, P700, and where $`|T=W_D^1_{j=3}^NT_{P700,j}|j`$ is a transient state representing the occupation probabilities right after a detrapping event has occurred. Substituting $`\mathrm{\Delta }`$, thus factored, into Eq. (39) the expansion terms for the average excitation lifetime can be written

$`\tau _1`$ $`=`$ $`q_\mathrm{𝟏}\tau _{soj},`$
$`\tau _2`$ $`=`$ $`q_\mathrm{𝟏}q_T\tau _{soj},`$
$`\tau _3`$ $`=`$ $`q_\mathrm{𝟏}q_T^2\tau _{soj},`$ (40)
$`\mathrm{}`$
where the defining parameters are

$`\tau _{soj}`$ $`=`$ $`\mathrm{𝟏}\left|\kappa ^1\right|T,`$
$`q_\mathrm{𝟏}`$ $`=`$ $`{\displaystyle \frac{W_D}{N1}}P700\left|\kappa ^1\right|\mathrm{𝟏},`$
$`q_T`$ $`=`$ $`W_DP700\left|\kappa ^1\right|T.`$ (41)
The terms appearing in Eq. (40) can be summed to give

$`\tau =\tau _0+{\displaystyle \frac{q_\mathrm{𝟏}}{1q_T}}\tau _{soj}.`$ (42)
The expansion decomposes the average excitation lifetime in terms of multiple escapes and returns (sojourns) to the reaction center, and hence, we call it the sojourn expansion.

Terms arising in the sojourn expansion for an effective Hamiltonian with full Coulomb couplings (see row 1 of Table 2) is provided in Table 3. It can be seen that the time for a first delivery of excitation to P700, approximately given by $`\tau _0`$, constitutes less than half of the average excitation lifetime at room temperature.

## 6 How optimal is the design of photosystem I?

At first glance the arrangement of chlorophylls in PSI appear to be random except for a pseudo-$`C_2`$ symmetry (see Fig. 1). Especially when compared to the regular circular aggregates found in purple bacterial photosynthetic systems, this brings up the question of whether the particular arrangement of chlorophylls in PSI is optimal in some way.

Even though the fitness landscape over which a photosynthetic apparatus should be judged is many dimensional, a crude measure of optimality is given by the quantum yield of the system. Ignoring all other concerns for optimality from an evolutionary perspective such as assembly, aggregation and interaction with other biomolecules in the cell, we examine whether the chlorophyll aggregate in PSI can be rearranged to result in a better quantum yield.

In order to keep the number of parameters to be varied at a manageable level, the chlorophyll positions (as defined by their Mg atoms) as well as the chlorophyll site energies are kept fixed. Instead we consider an ensemble of alternate PSI’s formed by random independent reorientations of constituent chlorophylls.

It is costly to repeat the semi-empirical computation of full Coulomb couplings between chlorophylls for a large ensemble of chlorophyll arrangements. Therefore, we revert to dipolar couplings for this study. The ensemble in question is formed by taking the vector for the transition dipole moment at each chlorophyll and multiplying it by a randomly generated $`SO\left(3\right)`$ rotation matrix, thus reorienting it. A new effective Hamiltonian is then computed using the couplings between these new transition dipole moment vectors as explained in Sec. 2, while keeping the same site energies as before. The effective Hamiltonians thus generated are used to compute the quantum yield and the average excitation lifetime at room temperature. For this study identical site energies (675 nm) for all chlorophylls are chosen except for P700, the site energies of which are chosen in such a way as to coincide with the observed absorption peak, as explained in Sec. 5. Therefore the results of this study should be compared with the original arrangement of chlorophylls given in row 2 of Table 2.

Fig. 11 contains a histogram of quantum yields and average excitation lifetimes for this ensemble at room temperature. The values for the original arrangement of chlorophylls are indicated by an arrow. It is seen that the original chlorophyll arrangement is nearly at the top of the distribution, though a few arrangements can be found to have a better quantum yield. It does not necessarily imply, however, that the efficiency of PSI can be easily improved upon. Not only is the study of quantum yields outlined above approximate in nature, but also there is no way to accurately model the constraints which had to be satisfied over the course of the evolution of this system. Nevertheless, it is impressive to find the original chlorophyll arrangement at near top efficiency.

The relatively small variation in the quantum yield across the ensemble seen in Fig. 11 brings up the question as to whether the apparent optimality is an artifact that will disappear in a more accurate computation. This small variation is a result of the nearly three orders of magnitude difference between the dissipation and charge transfer rates. As an ultimate test, full Coulomb couplings need to be computed for each system in an ensemble of randomly oriented chlorophyll aggregates and the effect of individual site energies also need to be taken into account. This is beyond the scope of this study due to prohibitively high computational costs. Nevertheless, as a consistency check, the ensemble described in Fig. 11 may be replaced by one with an alternate (random) set of site energies. It is seen that at room temperature the result of near optimal efficiency changes little under an assumption of different site energies (not shown). This is because the transfer rates at room temperature are influenced more strongly by couplings than by site energies, unless site energy variations are large enough to be comparable to the spectral width of individual chlorophylls (see Eq. 14). Also, the similarity between the networks of excitation transfer pathways for dipolar and full Coulomb couplings seen in Fig. 6 seems to suggest that the general properties of excitation transfer dynamics and therefore also the optimality will not be greatly effected by going from the dipolar to the full Coulomb couplings.

## 7 Conclusion

The availability of a high resolution structure of the cyanobacterial photosystem I provides new opportunities to investigate the mechanism of light harvesting in this large pigment-protein complex. In this study, the excitation transfer dynamics at room temperature has been examined starting from an effective Hamiltonian based on interactions between chlorophylls.

As a first approximation the effect of chlorophyll-protein interactions on the site energies has been ignored. However, we have used the available spectral data to construct a more refined Hamiltonian with limited information on site energies. Also, a prediction on the overall heterogeneity of site energies is made.

The calculated quantum yield and average excitation lifetime at room temperature compare favorably with experimental results. The sojourn expansion method resolves details of excitation migration in PSI, namely repeated approaches to and escapes from the reaction center. Studying an ensemble of Hamiltonians corresponding to copies of PSI with varying site energies, it is seen that the yield and the lifetime are influenced largely by inter-chlorophyll couplings and not so much by fluctuations of the site energies at room temperature. The sojourn expansion reveals that the time for a first delivery of excitation to P700 contributes less than half of the average excitation lifetime. At room temperature, about two thirds of the lifetime results from detrapping events.

However, the low temperature behavior of PSI is more difficult to describe accurately than the room temperature behaviour. Not only is the low temperature absorption spectrum impossible to reproduce without a knowledge of chlorophyll site energies, but also overlap integrals contributing to the excitation transfer rates become much more sensitive to the difference between the donor and the acceptor site energies due to the smaller linewidths at low temperatures. Furthermore, the broadening of lineshapes due to electron-phonon couplings also need to be carefully taken into account at low temperatures as this effect may have a significant impact on the overlap integrals in the framework of Förster theory. Ignoring these contributions results in an unrealistic picture of excitation transfer dynamics at low temperature.

The robustness and optimality of the chlorophyll arrangement in PSI has been examined by means of various computational experiments. On the one hand, it is seen that the efficiency of the chlorophyll network of PSI is robust against perturbations such as the pruning of individual chlorophylls or fluctuations of chlorophyll site energies. On the other hand, a study of an ensemble of alternate PSI’s with randomly reoriented chlorophyll aggregates shows that the original arrangement actually has a near optimal quantum yield.

## Note added in revision

After submission of this paper we have been informed of a recent calculation of chlorophyll site energies for cyanobacterial PSI . When these site energy values are substituted in the computations outlined above we obtain at room temperature an average excitation lifetime of 33.9 ps and a quantum yield of 96.6 $`\%`$, which compare favorably with the results given in this paper. The overall width of the distribution of site energies, including the red chlorophylls and P700, is around 260 cm-1, which is larger than the estimate made in this paper of 180 cm-1 based on the bulk of the low temperature spectrum. These new results support our suggestion that the room temperature excitation transfer dynamics of PSI can be described well without detailed knowledge of site energies. Repeating the optimality study discussed in section 6 also for the set of site energies in Ref.  yields results similar to those shown in Fig. 11.

## Acknowledgements

The authors would like to thank A. Damjanović, J. Gullingsrud, S. Hayashi, R. Knox, G. Small, E. Tajkhorshid, and R. van Grondelle for useful discussions. This work has been supported by NIH grants PHS 5 P41 RR05969 and PHS 1 R01 GM60946 and the Fetzer Institute.

## Appendix A Full Coulomb computation of inter-chlorophyll couplings

Here we outline the method employed to compute the electronic couplings, $`W_{\alpha \beta }`$, between two chlorophylls semi-empirically. This section follows Ref.  to which the reader is referred for a more detailed treatment.

The electronic coupling between a donor and an acceptor molecule can be written

$$W_{\alpha \beta }=\underset{i,jI_D}{}\underset{R,SI_A}{}C_{ij,RS}^c\times \psi _\alpha ^{}|^{00}O_j^i|\psi _\alpha \times \psi _\beta |^{00}O_S^R|\psi _\beta ^{},$$
(43)
where $`I_D`$ and $`I_A`$ denote the set of atomic orbital indices of the donor and acceptor chlorophylls, and $`C_{ij,RS}^c`$ describe the Coulomb integrals involving atomic orbitals labeled by $`i,j,R`$ and $`S`$. The spin tensors $`{}_{}{}^{00}O_{j}^{i}`$ and $`{}_{}{}^{00}O_{S}^{R}`$ prompt the intramolecular transitions $`|\psi _\alpha |\psi _\alpha ^{}`$, $`|\psi _\beta ^{}|\psi _\beta `$. The Coulomb interaction has zero rank, as denoted by the $`00`$ superscript. Therefore it only proceeds through singlet-singlet transitions. The Coulomb integral $`C_{ij,RS}^c`$ can be approximated  as $`S_{ij}\frac{e^2}{R_{ij,RS}}S_{RS}`$, where $`S_{ij}`$ and $`S_{RS}`$ are the atomic-orbital overlap integrals, and $`R_{ij,RS}`$ is the distance between the midpoint of atoms $`i`$ and $`j`$ and the midpoint of the atoms $`R`$ and $`S`$. As suggested in , $`S_{ij}`$ is taken to be $`1`$ when $`i=j`$, 0.27 when atoms $`i`$ and $`j`$ are joined by a chemical bond, and zero otherwise.

To evaluate the transition density matrix elements, one requires the description of the chlorophyll electronic states involved in the excitation transfer process. We choose a semi-empirical description for the chlorophyll electronic states as provided by the Pariser-Parr-Pople (PPP) Hamiltonian , which has been used earlier in . The PPP Hamiltonian

$`H_{PPP}`$ $`=`$ $`{\displaystyle \underset{i<j}{}}Z_iZ_jR_{ij}+{\displaystyle \underset{i,\sigma }{}}\left(I_i{\displaystyle \underset{ji}{}}Z_jR_{ij}\right)n_{i\sigma }`$ (44)
$`+`$ $`{\displaystyle \underset{ij,\sigma }{}}t_{ij}c_{i\sigma }^+c_{j\sigma }+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j,\sigma ,\stackrel{~}{\sigma }}{}}R_{ij}n_{i\sigma }n_{i\stackrel{~}{\sigma }}`$
involves orbitals of $`\pi `$-type only. The creation and annihilation operators, $`c_{i\sigma }^+`$ and $`c_{j\sigma }`$, act on mutually orthogonal atomic $`\pi `$-orbitals. The corresponding number operator is given by $`n_{i\sigma }=c_{i\sigma }^+c_{i\sigma }`$. $`R_{ij}`$ is the effective electron-electron repulsion integral between an electron in atomic obital at site $`i`$ and one in orbital at site $`j`$ and $`t_{ij}`$ denotes the resonance integral between atoms $`i`$ and $`j`$. $`I_i`$ is the effective ionization potential of an orbital at site $`i`$. $`Z_i`$ is the net charge of the core at atom $`i`$ which was chosen to be 1. The resonance integral, $`t_{\alpha \beta }`$, is evaluated according to the empirical formula

$$t_{\alpha \beta }=\gamma _0+3.21\left(r_{ij}1.397\mathrm{\AA }\right),$$
(45)
where $`\gamma _0`$ is a constant and $`r_{ij}`$ is the distance between the nuclear sites $`i`$ and $`j`$. The effective electron-electron repulsion integral, $`R_{ij}`$, is calculated through the Ohno formula ,

$$R_{ij}=14.397eV\left[\left(\frac{2\times 14.397eV}{R_{ii}+R_{jj}}\right)^2+\frac{r_{ij}^2}{\mathrm{\AA }^2}\right]^{\frac{1}{2}}.$$
(46)
Assuming the same values of the semiempirical parameters as in Ref.  (listed in Table 4), we perform a SCF-CI calculation, including all single excited $`\pi `$-orbitals, for the singlet states of a chlorophyll.

The chlorophyll electronic structure calculations are based on geometries of a simplified chlorophyll analog, displayed in Fig. 12 (solid line). In this analog, the double bond of ring II (dotted line in Fig. 12) is taken out. Thus the chlorophyll is almost symmetric about the magnesium atom except that there is a carbon atom at one end and an oxygen atom at the opposite end. For the calculations performed on the analog structure, the $`Q_x`$ and $`Q_y`$ states are easy to identify. However, for the calculations based on the real chlorophyll structure, where the electronic states are mixed, this identification is difficult. A fully symmetric analog, with both oxygen atoms at the end, was also proposed in  based on the same consideration and the error arising from the simplified chlorophyll analog was shown to be insignificant compared to the systematic errors.


text:
# More on the determination of the coronal heating function from Yohkoh data

## 1 Introduction

The most recent solar X-ray and UV missions (Yohkoh, SoHO, TRACE) have produced a significant amount of spectacular data and images of the corona. The improved spatial resolution and sensitivity of the detectors have revealed the highly structured and variable nature of the corona and have stimulated the interest of solar physics community.

Whereas the morphological interpretation of the images is often straightforward and has highlighted the presence of a variety of new structures and phenomena, the physical interpretation of the data and the inference of indications about items such as the heating of the corona is a much more delicate matter. A possible problem is that many physical and instrumental effects act concurrently and non-linearly. The coronal plasma is highly thermally conductive along the magnetic field lines but much less across them. As a result, non-uniform and localized heating sources which ignite coronal loops will at the same time be elusive, because the temperature increases rapidly along the whole loop, and determine a significant thermal stratification across the loop and along the line of sight. Since the coronal plasma is mostly optically thin, a telescope observing a loop will then collect the emission from multiple thermal components along the line of sight, filter them with different weights, depending on its passband, and sum them up to yield a single data number for each image pixel. The reconstruction of the original thermal structure of the plasma column in a pixel is a difficult task and some attempts have been made possible by multi-line observations with SoHO (Schmelz et al. 2001). When analyzing data from wide-band imaging instruments such as Yohkoh/SXT, one approach is to obtain maps of the weighted average temperature of the various thermal components along the line of sight. The interpretation of temperature distributions is somewhat easier along bright loop structures, because there one has the reasonable expectation that one thermal plasma component has a much larger emission measure than the others, and therefore dominates the emission and can be determined reliably.

This vision has driven recent studies of the thermal structure along coronal loops, and of its usage as indicator of the distribution of the heating deposition within the loop (Priest et al. 2000, Aschwanden 2001). The concept of these works was to compare the temperature distributions of bright coronal loops derived from the data with those obtained from models of thermally structured single coronal loops (or combinations of them). Priest et al. (2000) used an analytical approach and devised a method to infer the heating distribution along the loops from the analysis of temperature profiles, showing that along several bright steady-state large-scale loops observed at the solar limb the temperature distribution indicates uniformly distributed heating. Aschwanden (2001) revisited the analysis of one of the loops selected by Priest et al. (2000), and found that the previous approach led to inconsistent results and, in particular, to an unreasonably large plasma column depth required to match the emission measure inferred from the data. Using a two component loop model to explain the observations, he also found that the data are consistently fitted instead by a model coronal loop with a heating highly localized at the footpoints and with a cool background component.

In this work we present a further analysis of the same loop system, motivated as follows. The selected loop system is relatively bright, well outstanding at the solar limb. We might then suppose relatively large emission measure, high plasma pressure, and assuming that it is a steady structure, we can also infer, according to loop scaling laws (Rosner et al. 1978), high temperature at their apex, maybe higher than $`2`$ MK, which is the temperature obtained from the previous data analyses. On the other hand, neither previous analyses include, or include only partially, a component which instead may be important: “the loop system is embedded in a hazy background” (Aschwanden 2001). This background is significant, especially in the region of the loop apex, where the loop signal is fainter than elsewhere, due to the plasma gravitational stratification. One possible hypothesis, that we pursue here, is that the haze may be the effect of the presence of many other disordered faint loops, or of a streamer, which extend over the whole loop system and intersect the analysed structure along the line of sight. Such background may not be described as a simple hydrostatic equilibrium, and deserves proper attention, since it may affect systematically the filter ratio and temperature values. Another item stimulating our analysis is the finding by Aschwanden (2001) of a solution with a very small scaling height of the heating and temperature maxima at the loop footpoints, quite a puzzling result in the light of previous well-established loop modeling (Rosner et al. 1978, Serio et al. 1981).

Here we revisit the loop system, describing the main loop component with conventional hydrostatic single loop models and deriving the background component directly from the data. We compare the brightness profiles observed along the loop in the two SXT filter passbands with and without background subtraction to those obtained from the loop models. The modeling explores three different scenarios of heating deposition: uniform distribution, localized at the footpoints and localized at the apex, i.e. a more complete sampling than in Aschwanden (2001). The model results are compared to data with detailed fitting procedures, and attention is paid to the values of column depth obtained to match the observed fluxes. Section 2 describes the analysis of the data and their comparison with loop models; Sect. 3 discusses the results and their implications.

## 2 Data Analysis and Modeling

### 2.1 The Data

The loop system was observed with Yohkoh/SXT on 3 October 1992. The data consist of a sequence of full-disk 512$`\times `$512 pixel images ($`5`$” pixel side), basically in the Al.1 and Al/Mg/Mn filter passbands, taken with various exposure times during the whole day. The data have been preliminarily processed with the standard Yohkoh data analysis system (SXT\_PREP).

We have checked that the loop to be analyzed remains quite stable and stationary during the whole day. In order to enhance the signal-to-noise ratio, we then average the long-exposure (2.7 s and 15 s) images over the whole day (as done in Aschwanden 2001, 57 images in the Al.1 filter and 55 images in the Al/Mg/Mn filter). Fig. More on the determination of the coronal heating function from Yohkoh data shows the resulting image in the Al.1 filter passband. We extract the loop flux profiles within the arch-like strip (14 pixels wide) outlined in the figure. The strip has been divided into 31 sectors with approximately equal area. The flux profiles are then the sequences of Data Numbers averaged over each sector along the loop. DN fluctuations from each image to the other have been taken as the main source of data uncertainty.

As remarked by Aschwanden (2001), the loop system is embedded in a ”hazy” background. Here we make the hypothesis that this background is mostly due to the presence of other disordered faint loops and/or of a streamer which extend over the whole loop system and intersect the loop along the line of sight. Since the brightness of the haze is comparable to the loop brightness in the upper loop region, the effect of the background may be significant in the analysis and modeling of the loop emission. Some enhancement of the data quality was achieved by Priest et al. (2000) by reducing the radiation scattering due to SXT wings (Foley 1998), but, as we have checked (using the procedure SXT\_SCATTER), the amount of scattering ($`20`$ %) cannot account for the high background.

In the hypothesis that the background emission includes disorganized contributions from other faint structures, our approach has been to evaluate it with a model-independent procedure and to subtract it from the emission measured along the loop. To evaluate the background, we measure the emission in a region below the main visible arcade, where the emission is relatively faint, unstructured and as far as possible free of small evident loops. We then assume that, in the loop region, the background is stratified in the radial direction from the Sun center, and extends to involve the whole loop system. We extract the background fluxes (one for each filter) along the rectilinear strip (ten pixels wide) outlined in Fig. More on the determination of the coronal heating function from Yohkoh data, roughly perpendicular to the disk limb. This strip is chosen to have a signal decreasing upwards as far as possible smoothly. The resulting radial background profiles are then subtracted from the respective loop flux profiles, as a function of the distance of the loop sectors from the limb. One of the background profiles is shown in Fig. More on the determination of the coronal heating function from Yohkoh data, described below. The background flux results to be negligible nowhere along the loop, being of the order of 50% of the total flux. This procedure also allows us to limit the number of free parameters in fitting the dominant loop structure, which is the main target of the analysis.

We extract the loop flux within the arch-like strip shown in Fig. More on the determination of the coronal heating function from Yohkoh data which follows the geometry of the loop system. The loop enclosed in the strip does not appear to be perfectly semicircular, and symmetric around an axis perpendicular to the limb. We have assumed this to be due to a projection effect: the loop may not lie on a plane perpendicular to the line of sight. Then we have considered the data along half of the loop (the right leg) and corrected the coordinates along the loop for the projection effect (assuming a rotation around the symmetry axis only). As shown by the modeling results below, this allows us to make detailed comparisons with semicircular symmetric loop models.

### 2.2 Modeling

As done in Priest et al. (2000) and Aschwanden (2001), we try to compare the loop data with hydrostatic loop models. Priest et al. (2000) compare temperature values derived from the data with those predicted by models. Temperature values are obtained from the data with the filter ratio method, i.e. using the ratio of the flux in two filter passbands. The temperature values are then affected by two sources of uncertainties (one for each filter Data Number) and by the presence of high background. Instead of using temperature, our approach is then to compare flux profiles synthesized from model results directly to observed flux profiles corrected for background subtraction. The diagnostics on plasma temperature will then descend from the parameters of the model (or models) which is (are) found to provide a good description of the data.

The loop system is stable and steady on a time scale ($`1`$ day) much longer than the relevant physical time scales, i.e. the sound crossing time and the cooling times ($``$ few hours). We then consider standard hydrostatic single loop models as in Serio et al. (1981), semicircular, with constant cross-section along the loop and symmetric with respect to the loop apex. Given the extension of the loop system, gravity cannot be ignored and the loop is assumed to be perpendicular to the solar surface. The loop half-length is set to L = 380 Mm as in Aschwanden (2001). The model is based on the equations of plasma hydrostatic equilibrium and of energy balance among radiative losses, a phenomenological term of local heat input in the plasma and thermal conduction according to Spitzer’s formulation. The lower boundary of the loop is set at $`2\times 10^4`$ K, and there we assign the base pressure and the conduction flux ($`F_C0`$).

We generate models with various loop base pressures, corresponding to loop apex temperature logarithmically spaced (with step $`0.05`$) in the range 1 – 10 MK, for three different distributions of the heating along the loop: uniform heating, heating localized at the footpoints, and heating localized at the loop apex. The heating at the footpoints is described with a heating distribution exponentially decreasing upwards from the loop base, with a scale height $`s_H=127`$ Mm, i.e. 1/3 the loop half length. With this heating distribution we explore solutions for loop maximum temperatures below $`5`$ MK. We do not consider models with significantly smaller heating scale heights: they lead to loops with temperature maximum localized in the region of the loop footpoints. Such loops have been shown to be unstable (Rosner et al. 1978, Serio et al. 1981), and therefore unlikely to describe the loop analyzed here, which appear to be steady for a time of the order of a whole day. For the cases of heating at the apex, we assume a Gaussian distribution centered at the apex, with width $`\sigma _H=10`$ Mm.

From the temperature and density profiles obtained from the models we synthesize the loop emission filtered through the two SXT filter passbands, with the same procedure as in Aschwanden (2001). We then perform a best-fit procedure to the data along half of the loop, by finding analytically the normalization factor which minimizes the $`\chi ^2`$ of each model to the data points. This normalization factor is the only parameter of the fitting and represents the plasma column depth of the loop system along the line of sight. We then look for the models yielding the lowest $`\chi ^2`$ values among the models with different maximum temperatures and the same heating distribution, and then for the best-fit model among those with different heating distributions. This procedure is repeated on the data, both with and without background subtraction. In order to make the fitting in some way independent of possible errors in the determination of the loop length and of the real position of the relevant footpoint (see also the problem of submerged footpoints presented in McKay et al. 2000), the first three data points from the left are excluded from the $`\chi ^2`$ calculation, and we repeat the fitting for two different values of the relative shift of the model to the data, one assuming the model loop footpoint coincident with the data origin, and the other shifted backwards by $`2\times 10^9`$ cm. In the following we will show results only for the value of the shift yielding the better $`\chi ^2`$ value. We cannot exclude also some effects on the results due to loop opening at the base the loop. These are not expected to be important in corona (Serio et al. 1981), and may somewhat influence the normalization factor (Aschwanden 2001), but warnings on this have been raised due to the possibility that the coronal footpoint of a loop may have any base temperature since it may virtually be thermally isolated from the chromosphere (Dowdy et al. 1985).

Table 1 and Figs. More on the determination of the coronal heating function from Yohkoh data and More on the determination of the coronal heating function from Yohkoh data show some relevant results of the modeling, separating the cases of data without (hereafter unsubtracted data) and with background subtraction, and of models with uniform, base and top heating. Table 1 shows the significant parameters of the models in Figs. More on the determination of the coronal heating function from Yohkoh data and More on the determination of the coronal heating function from Yohkoh data, and in particular the pressure at the base of the model loop, its maximum temperature (at the top of the model loop), the (maximum) heating rate per unit volume, the loop plasma column depth derived as normalization factor, the reduced $`\chi ^2`$ values (11 d.o.f.), respectively obtained for the two filters, and the shift value yielding the best $`\chi ^2`$.

In Fig. More on the determination of the coronal heating function from Yohkoh data, each panel shows the data (data points) in both filters, the best-fitting model (solid lines) and other two models, one with lower and one with higher maximum temperature, for comparison. The origin of the reference system is at the loop footpoint. The panels on the left show data without background subtraction, those on the right data with background subtraction. From top to bottom, Fig. More on the determination of the coronal heating function from Yohkoh data shows fitting results for models with uniform heating, heating at the footpoints and heating at the apex, respectively. The upper left panel also shows one of the background profiles which are subtracted from the data to obtain the values shown in the lower two panels. The (similar) profile in the other filter passband is not shown for the sake of clarity.

Fig. More on the determination of the coronal heating function from Yohkoh data shows plots of the minimum $`\chi ^2`$ obtained from fitting the data with models with different maximum temperature for the two filters (data and models order as in Fig. More on the determination of the coronal heating function from Yohkoh data).

Several considerations are in order. The trends and values of the unsubtracted data are similar to those shown in Aschwanden (2001). Due to their lower statistics, the error bars of the background-subtracted data are larger than those of the unsubtracted data, and they are largest close to the loop apex, and in the Al/Mg/Mn filter passband. The background profile is monotonic, decreases toward the loop apex, and its trend is only partially similar to that of the full loop signal along the loop. The background-subtracted flux becomes quite small around the loop apex, below 1 DN/s in the Al/Mg/Mn filter passband.

Loop models with uniform heating resemble the trends of the unsubtracted data more closely than those of the background-subtracted data, although the latter ones yield better $`\chi ^2`$ values (because of the lower data statistics). The loop model with uniform heating best-fitting both the unsubtracted data and the background-subtracted data is the one with maximum temperature 3.3 MK (model 2 and 11 in Table 1) and a pressure $`0.4`$ dyne cm-2. The values of column depth obtained in the two filters are similar and reasonable ($`20`$ pixels), but the fitting is poor (high $`\chi ^2`$) in the region around the loop apex (Fig. More on the determination of the coronal heating function from Yohkoh data). Notice that the 2 MK models (e.g. model 1 and 10 in Table 1) yield an even worse fitting in the same region, and too large column depths, in agreement with Aschwanden (2001).

Models with heating at the footpoints are more in agreement with data both with and without background subtraction. Unsubtracted data are best-fit by the model with temperature 2.2 MK at the top (base pressure 0.2 dyne cm-2, model 5 in Table 1). However, the $`\chi ^2`$ value is high, and the column depth is unreasonably large (larger than the solar radius). Background-subtracted data are best-fitted by a model with a lower maximum temperature (1.9 MK, model 14 in Table 1), which yields a better (but still not statistically acceptable) $`\chi ^2`$ value, but an even larger column depth. Notice that both best-fit base-heated models, as well as those with uniform heating, do not describe well the flux at the footpoint, especially in the softer Al.1 filter passband, and systematically overestimate the flux at the loop apex.

Going to models with heating at the loop top, the one best-fitting the unsubtracted data (2.9 MK, model 8 in Table 1) largely underestimates the emission in the apex region, and yields very high $`\chi ^2`$ values and very large column depth values. On the other hand, the background-subtracted data appear to be adequately described in both filters by the (best-fit) top-heated loop model with a maximum temperature of 3.7 MK (model 17 in Table 1): the observed flux is well described anywhere along the loop, at the footpoint and at the apex altogether; the values $`\chi _\nu ^21`$ make the fitting statistically acceptable in both filters, and make such model the absolute best-fitting one among all the models explored here. Furthermore, the column depth values ($`10^{10}`$ cm) are reasonable for large scale or arcade structures. The low emission obtained at the apex with the best-fit top-heated model is the natural result of the localized heating excess, which steepens the temperature distribution to be more peaked at the apex, and makes the density decrease to maintain the pressure balance (Fig. More on the determination of the coronal heating function from Yohkoh data). The maximum temperature and the pressure at the footpoint (0.3 dyne cm-2) are not unreasonably high for large scale loops.

The fitting of the unsubtracted data invariably shows well-defined $`\chi ^2`$ minima (Fig. More on the determination of the coronal heating function from Yohkoh data, left column). This is true also for background-subtracted data, with the exception of the fitting with models with uniform heating, which shows quite a constant $`\chi ^2`$ with model maximum temperatures, and similar $`\chi ^2`$ values at low ($`2`$ MK) and high ($`>5`$ MK) temperatures (but with very different column depth values). High temperature top-heated models fitting the background-subtracted data is the only combination to yield statistically acceptable results.

Notice in Fig. More on the determination of the coronal heating function from Yohkoh data that fitting unsubtracted data in the Al/Mg/Mn filter passband yields generally better results than in the Al.1 filter passband. The opposite is true for background-subtracted data.

## 3 Discussion and Conclusions

This work has been stimulated by the contrasting results obtained from the analysis of the same loop system observed with Yohkoh/SXT by two subsequent works. In the earlier one (Priest et al. 2000) the data, and in particular the temperature distribution along the loop, were best matched by a loop model with a uniform heating distribution. The later one (Aschwanden 2001) found instead that a model dominated by a heating deposited at the loop footpoints, with a minor cooler component, provides more self-consistent results, because it predicts also reasonable values of the plasma column depth.

The present work extends the analysis of the same loop system, by including models with heating at the loop apex, and by considering a significant background emission possibly due to the intersection with other structures, such as disorganized fainter loops and/or a streamer, along the line of sight, which adds to the loop emission. The analysis shows that, if background emission is not subtracted, none of single loop models explored acceptably fits the data, indicating that the system cannot be described only as in a simple hydrostatic equilibrium. A base-heated loop model with maximum temperature 2.2 MK yields the best results on unsubtracted data, but, besides the high $`\chi ^2`$, it involves a column depth larger than the solar radius and trends not matching the data in some regions along the loop. We did not explore loop models with heating confined in a narrow region near the loop footpoints, as in Aschwanden (2001): such heating leads to a temperature maximum localized at the loop footpoints, and such configuration is well-known to be unstable (Rosner et al. 1978, Serio et al. 1981). On the other hand, the coronal plasma density of $`10^{10.3}`$ cm-3 predicted by the base-heated “best” model loop in Aschwanden (2001) is unusually large for large scale coronal structures (e.g. Vaiana & Rosner 1978). The best-fit uniformly-heated loop (with maximum temperature 3.3 MK) matches the data slightly worse than the best-fit base-heated model, but implies much more reasonable values of the column depth.

We do find a model fitting acceptably the background-subtracted data: a loop heated at the loop apex, with a maximum temperature of 3.7 MK (model 17 in Table 1). The corresponding column depth ($`1/10`$ of the loop length) represents a reasonable aspect for a loop system. Although the lower statistics of background-subtracted data should allow in general to obtain more easily an acceptable fitting, the model above is the only one to fit the data with high statistical significance. The combination of an acceptable fitting and of a realistic column depth provides a scenario globally self-consistent and consistent with the data.

Several considerations are now in order:

* Our best model is a loop heated at the apex, with a relatively high (still reasonable) temperature, a pressure of 0.3 dyne cm-2 at the footpoints and a density decreasing to $`3\times 10^7`$ cm-3 close to the apex. Such pressure and density values are not unusual for large scale structures (e.g. Vaiana & Rosner 1978). A loop heating localized in corona is not a new result (e.g. Reale et al. 2000).
* The best fitting results have been obtained after subtraction of an unstructured background emission, extrapolated directly from the data, in the hypothesis that it is due mostly to chance alignment of disordered faint loops and/or of a streamer along the line of sight, and therefore not to be modelled in hydrostatic equilibrium. This is one possible description of the scenario, not necessarily the best one, but, in our opinion, it is at least as realistic as other more local or model-dependent estimations of the background.
* High data statistics, obtained in this case by averaging over tens of exposures, is important to discriminate among loop models with different temperature and heating location. Fig. More on the determination of the coronal heating function from Yohkoh data shows that larger error bars would made such discrimination much harder. This shows, once again, that the task of fitting data with loop models is not at all trivial, and requires both detailed modeling and high quality data.

Of course, this work suffers from limitations. The background evaluation may not be unique, and a single loop model may be a simplified approximation to describe part of an arcade. Also the exploration of the loop parameter space is forcedly limited, and we cannot exclude that a more refined tuning of the free parameters, and including other effects, such as strong geometrical variations at the base of the loop, could bring to equally good or better fitting results. However, our analysis and modeling lead to a statistically acceptable description of the data, and to physically sound results (loop heating and plasma conditions, column depth), and appears therefore to be adequate to the data quality.

This work shows that the analyzed observational data can be described with a conventional hydrostatic and non-isothermal loop models and with an unstructured vertically stratified spurious emission. The deposition of heat at the apex of large scale structures is a different result from those of both previous works, and may not be in agreement with recent results based on the analysis of data from the TRACE mission (Aschwanden et al. 2001). Since there are serious independent arguments which make the analysis and temperature diagnostics with narrow band EUV instruments a very delicate issue (Schmelz et al. 2001), we believe that, at variance with Aschwanden (2001), the problems of the diagnostics of the heating function in corona and of the coherent interpretation of observations with different telescopes (e.g. Yohkoh and TRACE) are both still open.

We do not pretend to put a conclusive word on this topic, nor claim that our results are necessarily the best ones, rather we want to stimulate the community toward a more careful attention to aspects of the data analysis, such as the evaluation of background emission, important for a correct interpretation and modeling of the data. Other authors (McKay et al. 2000) come to similar conclusions, although with a totally different analysis. The present work also emphasizes the importance of the selection of the loops to be analyzed: loops embedded in crowded and complex coronal regions may not be the best ones, because emission from structures of different kind and characteristics may intersect along the line of sight, which may be difficult to evaluate. Further extensive and systematic analysis and modeling of bright, and as far as possible isolated, loop structures may provide more constraints on the issue of coronal heating deposition.

I would like to thank Salvatore Orlando, Paola Testa and Salvatore Serio for their helpful criticism and interest in this work. This work was supported in part by Agenzia Spaziale Italiana and by Ministero dell’Istruzione, dell’Università e della Ricerca.