text:
# Origin of Scaling Behavior of Protein Packing Density: A Sequential Monte Carlo Study of Compact Long Chain Polymers

## I Introduction

Geometric considerations have lead to important insights about protein structures Chothia81 ; Chothia90 ; Maritan00\_Nature ; Banavar02\_Proteins ; Bagci02\_JCP . Voids are simple geometric features that represent packing defects inside protein structures. For multisubunit proteins such as GroEL and potassium channel, voids or tunnels of large size are formed by the spatial arrangement of multiple subunits, and are essential for the biological functions of these proteins Sigler98\_ARB ; Doyle98\_Science . In this study, we focus on voids formed due to packing defects that are not directly involved in protein function. For this purpose, we choose to study only structures of single domain proteins. Although these proteins are well known to be compact Richards94\_QRB , and their interior is frequently thought to be solid-like Hermans61\_JACS ; Richards97\_CMLS , recent calculations showed that there are also numerous voids buried in the protein interior LiangDill01\_BJ . The importance of tight packing in single chain protein is widely appreciated: packing is thought to be important for protein stability Eriksson92\_Science ; Privalov96 ; Dahiyat97\_Science , for kinetic nucleation of protein folding Ptitsyn98\_JMB ; PtitsynTing99\_JMB , and for successful design of novel proteins following a predefined backbone Dahiyat97\_Science . The conservation of amino acid residues during evolution may also be correlated with tightly packed sites Ptitsyn98\_JMB ; PtitsynTing99\_JMB ; MirnyShahk01\_ARBBS . In contrast, the potential roles of voids in affecting protein stability and in influencing tolerance to mutations and designability of proteins Li96\_Science ; Melin99\_JCP  are not well understood.

An important parameter describing packing is the packing density $`p_d`$, which is a quantitative measure of the voids and was first introduced to study proteins by structural biologists. This concept has been widely used in protein chemistry Richards94\_QRB ; Privalov96 . The scaling relationship of $`p_d`$ and chain length $`N`$ was first studied in reference LiangDill01\_BJ . $`p_d`$ can be thought of as the physical volume $`v_{vdw}`$ occupied by the union of van der wall atoms, divided by the volume of an envelope $`v_{env}`$ that tightly wraps around the body of atoms: $`p_dv_{vdw}/v_{env}`$ LiangDill01\_BJ . Voids contained within the molecule will not be part of the van der Waals volume $`v_{vdw}`$, but will be included in $`v_{env}`$. Using geometric algorithms, $`v_{vdw},v_{env}`$ and $`p_d`$ can be readily computed for protein structures in the Protein Data Bank Liang98a\_Proteins ; Liang98b\_Proteins .

In this work, we further study the scaling behavior of packing density $`p_d`$ with chain length of single domain proteins and explore the determinants of the observed scaling behavior. We seek to answer the following questions: Is the scaling behavior of $`p_d`$ unique to proteins? Are proteins optimized during evolution to eliminate packing voids? We introduce two new packing parameters $`n_\alpha `$ (the alpha contact number) and $`z_\alpha `$ (the alpha coordination number). We show that $`n_\alpha `$ characterizes protein packing very well with a linear scaling relationship with the chain length, and that a widely used parameter, the radius of gyration $`R_g`$, characterizes protein packing poorly. To overcome the attrition problem of low success rate in generating compact long chain polymers, we develop an algorithm based on sequential Monte Carlo importance sampling and succeed in obtaining thousands of very compact long chain off-lattice polymers up to $`N=2,000`$. We demonstrate that the scaling behavior of $`p_d`$ for proteins can be qualitatively reproduced by randomly generated polymers with rudimentary constraints of $`n_\alpha `$. Our simulation studies lead us to conclude that proteins are not optimized to eliminate voids during evolution. Rather, voids in proteins are a generic feature of random polymers with a “reasonable” (as measure by $`z_\alpha `$) compactness.

The paper is organized as follows. In the Method section, we first describe briefly how $`p_d`$ and $`n_\alpha `$ are computed from the dual simplicial complex of protein structure, and introduce an off-lattice discrete model for generating random polymer conformations. We next describe the sequential Monte Carlo importance sampling and resampling techniques that allow us to generate adequate samples satisfying various criteria of $`n_\alpha `$. In the Results section, we begin with the characterization of void properties of proteins by both $`p_d`$ and $`R_g`$. We then show the linear scaling behavior of $`n_\alpha `$ found in proteins. The scaling behavior of $`p_d`$ of random polymers generated by sequential Monte Carlo with chain length is discussed later. We conclude with summary and discussion of our results.

## II Methods

Protein Data. To avoid complications of multichain and multidomain proteins, we examine the packing density of proteins of single domain proteins. We collect proteins from the Pdbselect database Hobohm94\_PS  that contains only one domain, as defined as single chains in the Scop database with one numerical label SCOP .

Dual Simplicial Complex, Alpha Coordination Number, and Packing Density. We use alpha shape to characterize the geometry of protein structure. Alpha shape has been successfully applied to study a number of problems in proteins, including void measurement, binding site characterization, protein packing, electrostatic calculations, and protein hydrations Edels95\_Hawaii ; Liang98a\_Proteins ; Liang98b\_Proteins ; Peters96\_JMB ; Liang98\_PS ; LiangDill01\_BJ ; Liang97\_BJ ; Liang98\_BJ . Briefly, we first obtain a Delaunay simplicial complex of the molecule from weighted Delaunay triangulation, which decomposes the convex hull of atom centers into tetrahedra (3-simplices), triangles (2-simplices), edges (1-simplices), and vertices (0-simplices). We then obtain the dual simplicial complex of the protein molecule by removing any tetrahedra, triangles, and edges whose corresponding Voronoi vertices, edges, and planar facets are not fully contained within the protein molecule Edels95\_DCG ; Edels98\_DAM . The edges between atoms that are not connected by bonds corresponds to nonbonded alpha contacts. The total sum of the number of edges for each atom is the total number of alpha contacts $`n_\alpha `$. It reflects the total number of atoms that are in physical nearest neighbor contact with other atoms. These atoms have volume overlap and their corresponding weighted Voronoi cells intersect. The alpha coordination number is $`z_\alpha n_\alpha /n`$, where $`n`$ is the total number of atoms in the molecule (see Figure 1). In our calculation, we only consider nonbonded alpha contacts. Details of the theory and computation of alpha shape and dual simplicial complex can be found elsewhere Edels94\_ACMTG ; Liang98a\_Proteins ; Liang98b\_Proteins ; LiangDill01\_BJ .

We follow previous work in reference LiangDill01\_BJ  and define packing density $`p_d`$ as:

$$p_d\frac{v_{vdw}}{v_{env}}=\frac{v_{vdw}}{v_{ms}+v_{voids}}$$

where $`v_{vdw},v_{ms}`$ and $`v_{voids}`$ are van der Waals volume, molecular surface volume and the void volume of the molecule, respectively Liang98a\_Proteins . Packing density is computed with a solvent probe radius 1.4 Å, as described in reference LiangDill01\_BJ .

Growth Model for Off-Lattice Random Polymers. We use a modified off-lattice discrete $`m`$-state model first developed in reference ParkLevitt95\_JMB  to generate self-avoiding walks (SAWs) in three dimensional space. All monomers are treated as balls with a radius of 1.7 Å. For monomers $`i`$ and $`j`$ that are not sequence near neighbors ($`|ij|>2`$), the Euclidean distance $`d(i,j)`$ between them must be greater than $`2\times 1.7`$ Å so they are self-avoiding. Sequence neighboring monomers are connected by a bond of length 1.5 Å.

We use a chain growth model to obtain conformation of polymer of specified length Rosenbluth55\_JCP . There are $`m=32`$ possible states where the next monomer can be placed. They are evenly distributed spatially on a sphere of radius 1.5 Å centered at the current monomer. We forbid the placement of the new monomer anywhere on a cap of the sphere with an angle $`<60`$  from the entering bond. This ensures that there are no unnatural acute sharp bond angles. The remaining sphere is divided into 4 strips, each may have different width but is of equal surface area. For the 32 possible states, we place uniformly 8 points at the midline of each of the 4 strips. Following Park & Levitt ParkLevitt95\_JMB , the coordinates of each state are parameterized by two angels $`\alpha `$ and $`\tau `$ for ease of computation. $`\alpha `$ is the bond angle formed by the $`i1,i`$ and $`i+1`$-th monomers. $`\tau `$ is the torsion angle formed by four consecutive monomers.

Approximately Maximum Compact Polymer. In addition, we generate polymers that are approximately maximum compact based on the face centered cubic (FCC) packing of balls of 1.7Å radii. For hard spheres, FCC packing has recently been proved to have the tightest packing FCC ; Sloane98\_Nature . Because the distance between two balls in canonical FCC packing is $`2\times 1.7=3.4`$ Å, which is greater than the bond length 1.5 Å, we shorten the distance along bonds connecting contacting balls of radius 1.7Å to 1.5 Å. This mimics the bond length of the model polymer. Unlike FCC packing of hard spheres, bonded monomers here are allowed to have volume overlaps. Additionally, there are some boundary effects because bonds connecting balls in different layer have a distance $`>1.5`$ Å. Although mathematically unproven, we conjecture that this artificially constructed polymer represents conformations of SAWs that have very close to maximum compactness.

The packing density of canonical FCC packing by our method is 0.74 LiangDill01\_BJ . As described earlier, although FCC packing contains no voids, there are packing crevices or dead spaces that do contribute to the calculation of $`p_d`$ by our definition LiangDill01\_BJ . In approximately maximum compact polymer, because the distance between bonded balls is shorter than that in FCC packing, $`p_d`$ can be as high as 0.80 for polymers with a range of chain length.

Importance Sampling with Sequential Monte Carlo. Since we are simulating compact conformations that resemble proteins, we need an efficient method to generate adequate number of conformations satisfying protein-like compactness criteria. Here we use a sequential Monte Carlo (SMC) chain growth strategy Liu&Chen98 ; LZC02\_JCP , which combines importance sampling and the growth method. The main steps are shown in Figure 3.

Denote the conformation of a polymer of length $`t`$ as $`(x_1,\mathrm{},x_t)`$, where $`x_i`$ is the 3-dimensional location of the $`i`$-th monomer. Starting with fixed initial location $`(x_1,x_2)`$, we grow polymers by sequentially adding one monomer $`x_{t+1}`$ to occupy one of the 32-states connecting to the last monomer $`x_t`$ of the current chain. The monomer $`x_{t+1}`$ is randomly placed according to a sampling probability $`g_{t+1}(x_{t+1}|x_1\mathrm{}x_t)`$. In this study, the following function $`g_{t+1}`$ is used. Let $`\omega `$ be one of the 32-states connected to $`x_t`$ that satisfies the self-avoiding criterion. First we intitialize the number of neighbors $`n_e(\omega )`$ to $`\omega `$ as 1, and the Euclidean distance from $`\omega `$ to the nearest neighbor monomer $`d(\omega )`$ to 6 Å. We then increment $`n_e(\omega )`$ by the number of existing monomers within a distance of 6.0 Å to $`\omega `$. Among these monomers, we identify the monomer $`x_s`$ that is the nearest neighbor with the shortest Euclidean distance $`d`$ to $`\omega `$. We require in addition that the sequence separation $`|st|>3`$ so $`x_s`$ and $`x_t`$ are not sequence near neighbors. The distance $`d(\omega )`$ is then replaced by the value of $`d`$. The sampling probability is set as:

$$g_{t+1}(x_{t+1}=\omega |x_1\mathrm{}x_t)e^{E^{}(\omega )/T^{}}$$

where $`E^{}(\omega )=\mathrm{ln}\frac{[d(\omega )]^c}{n_e(\omega )}`$ is an artificial “packing energy” favoring more compact conformations, and $`T^{}`$ is a pseudo-temperature controlling the behavior of sampling. Using this energy function, growth to position $`\omega `$ with close nearest neighbor (small $`d(\omega )`$) and a large number of neighbors within a 6 Å distance (large $`n_e(\omega )`$) is favored. Here the adjustable parameter $`c`$ is used to balance the effect of $`d(\omega )`$ and $`n_e(\omega )`$. $`T^{}`$ controls the importance of compactness. At low $`T^{}`$, conformations generated are compact, but at high $`T^{}`$, the compactness criterion becomes less important.

According to the sequential Monte Carlo framework, the importance weight $`w_{t+1}`$ for the sampled conformation $`(x_1,\mathrm{},x_{t+1})`$ is updated as:

$$w_{t+1}=w_t\frac{\pi _{t+1}(x_1\mathrm{}x_{t+1})}{\pi _t(x_1\mathrm{}x_t)g_{t+1}(x_{t+1}|x_1\mathrm{}x_t)}$$

where $`\pi _{t+1}(x_1\mathrm{}x_{t+1})`$ is the target distribution at $`t+1`$. With a set of weighted samples $`\{(x_1^{(j)},\mathrm{},x_n^{(j)}),w_n^{(j)}\}_{j=1}^m`$, statistical inference on the target distribution $`\pi _n(x_1,\mathrm{},x_n)`$ can be made using

$$E_{\pi _n}[h(x_1,\mathrm{},x_n)]=\frac{_{j=1}^mw_n^{(j)}h(x_1^{(j)},\mathrm{},x_n^{(j)})}{_{j=1}^mw_n^{(j)}}$$
(1)
for most of proper function $`h`$.

The Target Distribution. We wish to generate random samples of polymer with different compactness criterion. This is achieved by using a target distribution $`\pi _n`$ which is uniform among all SAWs satisfying a compactness constraint. The constraint is set as follows. First, for each chosen pair values of $`(T^{},c)`$, we use the function $`e^{E^{}(c)/T^{}}`$ to generate 500 random conformations as a trial run. Ignoring the importance weights, we calculate the mean alpha coordination number $`z_\alpha ^{}(n,T^{},c)`$ of all the generated conformations. Then we set the target distribution $`\pi _n^{}`$ as the uniform distribution of all SAWs satisfying $`z_\alpha (0.8z_\alpha ^{}(n,T^{},c),1.2z_\alpha ^{}(n,T^{},c))`$.

We then rerun a large simulation with the same $`(T^{},c)`$ parameters and harvest the conformations, using uniform distribution with no restriction on the intermediate target distribution $`\pi _t`$ but take the truncated distribution $`\pi _n^{}`$ as the final target distribution. The truncation is archived by discarding all generated conformations that does not satisfy the constraint. Typically, the truncation rate is very small ($`<0.1\%`$). The bias in sampling is fully compensated by proper weighting.

Resampling. Because it is easy to have self-avoiding walks to grow into a dead-end, we use resampling to replace dead samples or samples with small weight to improve sampling efficiency LZC02\_JCP . Intuitively, we check regularly during the chain growth process whether a particular chain is stuck in a dead-end, or is too extended, or has too little weight. If so, this chain is replaced by the replicate of another chain that has the desired compactness. Both duplicate chains will then continue to grow, and the final two surviving chains will be correlated up to the duplication event. Conformations of the monomers added after the duplication will be uncorrelated. This resampling technique targets our simulation to specified configuration space without introducing too much bias where conformations all have desired compactness (see Figure 4).

Although we found that the total contact number $`n_\alpha `$ is an excellent parameter for characterizing protein, its calculation involves expensive computation of weighted Delaunay triangulation and alpha shape. We decide to use $`R_g`$ as a surrogate parameter during resampling. For resampling, we use the empirical relationship $`R_g(n)=2.2n^{0.38}`$, where $`n`$ is the number of monomers in the polymer, as described in reference Skolnick97\_JMB . This relationship has been used as a constraint in NMR protein structure determination HuangPowers01\_JACS . We have the following pseudo code for resampling LZC02\_JCP :

| Procedure Resampling ($`m,d_s,R_t`$) |
| --- |
| // $`m`$: Monte Carlo sample size, $`d_s`$: steps of looking-back. |
| // $`R_t`$: targeting $`R_g`$. |
| $`k`$ number of dead conformations. |
| Divide $`mk`$ samples randomly into $`k`$ groups. |
| for  group $`i=1`$ to $`k`$ |
|  | Find conformations not picked in previous $`d_s`$ steps. |
|  |  | //Pick the best conformation $`P_j`$ |
|  |  | $`P_j`$ polymer with $`\mathrm{min}|R_gR_t|`$ |
|  | Replace one of $`k`$ dead conformations with $`P_j`$ |
|  | Assign both copies of $`P_j`$ half its original weight. |
| endfor |

Here $`d_s`$ is used to maintain higher diversity for resampled conformations. That is, conformation that has been picked in the past $`d_s`$ steps are not available for resampling.

After resampling, the samples with their adjusted weights remain to be properly weighted with respect to the original target distribution. We can then calculate the expected alpha coordination number $`z_\alpha `$, expected packing density $`p_d`$, or expected value of any other function $`h`$ using Equation (1). With these sampling and re-sampling strategies, we can successfully grow thousands of self-avoiding walks of chain length up to $`2,000`$ using a Linux cluster of 40 CPUs.

## III Results

Packing Density. Figure 5a shows the correlation of packing density $`p_d`$ with the number of residues $`N`$ in real proteins. Similar relationship has been observed in reference LiangDill01\_BJ . Here we further restrict the samples to be of single domain by Scop annotation SCOP . We found that $`p_d`$ decreases with chain length. That is, short chain proteins have high packing density $`p_d`$, but $`p_d`$ decreases from $`>0.85`$ to about $`0.740.75`$ when the chain length reaches about 190 residues. After reaching this length, proteins seem to be indifferent about the existence of voids. This suggests that maintaining high packing density is only characteristic of short chain proteins.

Radius of Gyration of Proteins. To identify the factors that dictate the scaling behavior of $`p_d`$ with residue number $`N`$, we need to determine whether such scaling is due to physical constraints of statistical mechanics or the product of extensive optimization by evolution. We study this problem by examining the scaling behavior of $`p_d`$ with $`N`$ in random chain polymers generated by computer.

Because of the enormity of conformational space, we focus on random polymers that resemble proteins in some rudimentary sense. One possible criterion is the radius of gyration $`R_g`$. This parameter has been widely used as a macroscopic description of protein packing. For single domain proteins, however, we found that there is substantial variance in $`R_g`$ for proteins of the same chain length (Figure 5b). Therefore, $`R_g`$ characterizes protein packing rather poorly, and is unsuitable as a criterion for generating protein-like polymers for our purpose.

Alpha Contacts. An alternative global description of protein structure is the total number of nonbonded alpha contacts $`n_\alpha `$ defined by the dual simplical complex of the protein. In Figure 6a we plot $`n_\alpha `$ against the total number of atoms in the molecule $`n`$. As discussed before, these contacts are identified by computing the dual simplicial complex of the molecule Liang98a\_Proteins ; LiangDill01\_BJ . The total number of contacts $`n_\alpha `$ scales linearly with $`n`$. It also scales linearly with the protein chain length (or residue number $`N`$, data not shown). Regression leads to a linear relationship of $`n_\alpha =4.28n432`$, with $`R^2=0.995`$. The alpha contact number $`n_\alpha `$ therefore provides a more accurate global characteristic of protein than radius of gyration $`R_g`$. This linear scaling relationship of packing related property is similar to other linear scaling relationships observed for protein, for example, of empirical solvation energy Chiche90\_PNAS , protein surface area and protein volume LiangDill01\_BJ  with chain length. It is interesting to note that the value of $`x`$-axis intercept for $`n`$ of the linear regression model suggests that the size of a minimum protein would be in the order of 100 atoms, or about 9–10 residues.

The details of the linear scaling relationship are further examined in Figure 6b. It is a replot of Figure 6a after normalization by $`n`$. It showed that for proteins with 1,000 atoms or more ($``$ 120 residues), the parameter alpha coordination number $`z_\alpha =n_\alpha /n`$ is a constant of about 4.2. For smaller proteins ($`n<1,000`$), $`z_\alpha `$ ranges from 2.5 to 4.0. A nonlinear curve fitting leads to the relationship $`z_\alpha =ab/n`$, where $`a=4.27\pm 0.03`$, and $`b=4.2\times 10^2\pm 26`$. We decide to use $`z_\alpha `$ as the criterion to select random polymers generated computationally for packing analysis.

Targeted Sampling of Random Chain Polymer. Figure 7 shows typical conformations generated with different $`(T^{},c)`$ parameters and the conformation of maximally compact polymer. Figure 8 shows the histogram of $`z_\alpha `$ of the conformations at length 2,000 without weight adjustment generated using different $`(T^{},c)`$ parameters. It can be seen that the histograms for different values of $`(T^{},c)`$ do not overlap. This feature demonstrated that with properly chosen $`(T^{},c)`$, we can efficiently generate random polymers with $`z_\alpha `$ within a targeted range.

Packing Density of Random Chain Polymer. Figure 9a shows the relationship of $`c_\alpha `$ associated with each pairs of $`(T^{},c)`$ as a function of chain length $`n`$. It also shows the $`c_\alpha `$ value for the maximally compact conformations. Figure 9b shows the relationship of $`z_\alpha `$ and $`n`$. Note that the targeted $`z_\alpha `$ generated with different $`(T^{},c)`$ parameters give rise to different $`z_\alpha n`$ scaling behavior. Because the coarse grained random polymers generated here lack side chains, they are fundamentally different from real proteins. We therefore have experimented with several $`(T^{},c)`$ value. We find that protein-like scaling can be obtained for a wide range of $`(T^{},c)`$ values.

Figure 9b shows the average packing density $`p_d`$ for all conformations satisfying the constraint specified by different $`(T^{},c)`$ values. Except maximally compact conformations, the scaling of $`p_dn`$ of all other sets of polymers is remarkably similar to that of protein (Figure 5a).

Conformations from Set 1 ($`(T^{},c)=(1.0,0.0)`$) are more extended, and have lower average $`z_\alpha `$ (e.g., Figure 7a). Because there are fewer voids, they also have high $`p_d`$. Conformations in Set 3 ($`(T^{},c)=(0.1,0.6)`$) are more compact and make more nonbonded contacts and hence have high $`z_\alpha `$ values (e.g., Figure 7c). They also form more voids, and therefore have lower $`p_d`$ values. Set 2 ($`(T^{},c)=(0.67,0.0)`$) are conformations whose properties are between those of set 1 and set 3 (e.g., Figure 7b).

The relationship between $`z_\alpha `$ and chain length $`n`$ can be characterized by a nonlinear equation $`z_\alpha =ab/n`$, similar to that of real proteins. The sets of $`a`$ and $`b`$ obtained by curve fitting are listed in Table 1. We emphasize that these randomly generated self-avoiding walks are fundamentally very different from proteins: all residues are of uniform size, there are no side chains, and there is no hydrophobic or any other type of physical interactins in these polymers. Because it is impossible to quantitatively define a similarity metric that measures how different these polymers are from proteins, we are not able to decide which specific values of $`(T^{},c)`$ are optimal for modeling protein packing. Nevertheless, the scaling of $`p_dn`$ for all $`(T^{},c)`$ values is qualitatively quite similar to that of real proteins.

The relationship between $`p_d`$ and $`z_\alpha `$ at chain length 1,800 for self avoiding walks generated with different parameters $`(T^{},c)`$ are shown in Figure 10. The $`p_d`$ of both extended and maximally compact conformations have high $`p_d`$ values, but conformations with an intermediate value of $`z_\alpha `$ contain voids and have smaller $`p_d`$ values. This is similar to the relationship of $`p_d`$ and a compactness parameter $`\rho `$ (equivalent to $`z_\alpha `$ used here) studied in two dimensional lattice (see Figure 8 in reference LZC02\_JCP ).

## IV Summary and Discussion

It is well acknowledged that protein has high packing density $`p_d`$, as high as that of crystalline solids Richards94\_QRB . However, recent study suggested that there are numerous voids and pockets in proteins LiangDill01\_BJ . It was also found that about 1/3 of the residues in a protein deviates from the FCC close packing and have random positions Bagci02\_JCP . The simulation results presented here indicate that chain connectivity, excluded volume, and global compactness are the main determinants of the scaling behavior of voids and chain length in proteins. Unlike maximally compact polymers which maintains high packing density at all chain length, proteins and simple near compact polymers have large $`p_d`$ values only for relatively short chains. When the chain length reaches 190 residues for protein and about 600–700 for chain polymers, proteins and polymers have lower packing density and are quite tolerant to the formation of voids.

The global compactness is a necessary condition for the observed protein-like scaling behavior. However, not all parameters related to voids and compactness are equally appropriate. The data shown in Fig 6 suggests that the alpha coordination number $`z_\alpha `$ reflects basic intrinsic compactness properties of protein, which is absent in the widely used parameter $`R_g`$, the radius of gyration. The advantage of parameters such as $`z_\alpha `$ emphasizes the importance of accurate description of protein geometry and structure.

The parameters $`p_d`$ is biased towards short chain proteins. By definition, a polymer formed by 2, 3 or a small number of monomers do not have long enough chains to form voids, therefore all will have $`p_d=1.0`$. When chain length becomes longer, voids appear. Similarly, $`z_\alpha `$ is also biased towards short chains. For very short chain polymers where the chain has few turns, few non-bonded contacts exist and no voids are formed. In this case, $`z_\alpha `$ is low and $`p_d`$ is high. However, this small size effect disappears rapidly for our model conformations for maximally compact polymers. Small size effect therefore does not fully account for the scaling behavior of $`p_d`$ and $`z_\alpha `$ in proteins and in simulated random polymers.

There are major differences between self-avoiding walks we generated and real protein structures. Our SAWs have no side chains, and belong to the coarse-grain model where one monomer is represented as a ball. In addition, the target distribution of SMC sampling is the truncated uniform distribution of all geometrically feasible conformations. The truncation required is that polymers must satisfy a prescribed $`z_\alpha n`$ relationship. No physical forces such as hydrophobic interactions is used in the target distribution.

In this paper, we describe a novel approach to overcome the attrition problem in generating long chain compact self-avoiding walk. With sequential importance sampling and resampling, we have developed an algorithm that effectively sample rare events, i.e., compact self-avoiding walks. Success in generating thousands of off-lattice self-avoiding walks satisfying various desired compactness requirement is essential for studying the scaling behavior of $`p_d`$ in random off-lattice self-avoiding walks.

The main result of this paper is that with sequential Monte Carlo techniques it is not difficult to reproduce protein-like scaling behavior of packing density $`p_d`$ and chain length $`n`$ in generic chain polymers. With the guidance of rudimentary requirement of $`z_\alpha `$, this can be achieved under a wide range of $`z_\alpha n`$ relationships. We therefore conclude that proteins retain the same packing property of generic compact chain polymers. We further conclude that proteins are unlikely to be optimized by evolution to eliminated packing voids. This is in support of the insightful comments of Richards who suggested that an appropriate level of under-packing would be important for evolution to occur through random mutations Richards97\_CMLS .

Our study showed the importance of generic geometric packing related to $`z_\alpha `$ in reproducing protein-like $`p_dn`$ scaling behavior. To test further the role of geometric packing, the next step would be to examine the $`p_dn`$ scaling by generating more realistic compact random polymers with perhaps monomers of different sizes to model the side chain effects. Furthermore, with sequential Monte Carlo and other advanced sampling methods, various models of explicit side chains can be attached to main chain monomers. In addition, one could introduce various alphabet sets for the residues (such as the HP model) and corresponding potential energy function $`H`$. In this case, the target distribution can be the Boltzmann distribution $`\pi \mathrm{exp}(H/T)`$ instead of the uniform distribution of all SAWs. It would also be interesting to examine the $`p_dn`$ scaling of polymers of random sequences and of protein-like sequences with low energy in compact states.

## V Acknowledgments

We thank Profs. Herbert Edelsbrunner and Luhua Lai for valuable discussions. This work is supported by funding from National Science Foundation CAREER DBI0133856, DBI0078270, DMS0073601, CCR9980599 and American Chemical Society/Petroleum Research Fund 35616-G7.


format: text

text:
# Polarisation Observations at Jodrell Bank Observatory

## 1. Introduction

Polarisation can prove to be an invaluable tool in determining, among other things, the geometrical properties of objects and phenomena. A large amount of information can be retrieved from studying pulsar polarisation profiles. It can also be used to probe the emission mechanism and geometry of pulsars, pulsar precession and the magnetic field of our Galaxy.

The Jodrell Bank database of polarisation observations covers more than 16 years and over 300 pulsars observed at frequencies of 230, 400, 600, 920, 1400 and 1600 MHz (Gould & Lyne 1998), providing an ideal archive for thorough polarisation research.

## 2. Necessity of calibration: the case of Jodrell Bank

Acquiring accurate polarisation profiles is not a trivial task and a considerable number of calibration issues needs to be addressed before the corrected profiles can be actually studied (e.g. Hamaker, Bregman, & Sault 1996 or Stinebring et al 1984). Some of the main problems are:

* Gain calibration
* Feeds’ non-orthogonality and ellipticity
* Parallactic angle correction due to Earth’s rotation
* Faraday depolarisation due to the ionosphere & the ISM
* Signal path length difference

The calibration process at Jodrell Bank is based on a noise diode. In particular, information from the 4 polarimeter channels with the diode being alternatively switched on and off is used. This is combined with observations of an unpolarised continuum source of known flux density (e.g. 3C123) and observations of 1 degree off the pulsar and at half power.

The profiles resulting from this procedure are gain calibrated, they have any channel offsets subtracted and they are corrected for parallactic angle rotation as well as the noise diode position angle. However, effects related to the feeds misalignment and depolarisation from Faraday Rotation and path differences in the signal may still produce intrumental errors.

Therefore, although the calibration already performed is adequate for timing studies, more calibration observations are required to correct for the above effects as well. Such observations are not always available in the archive though and in order to minimise the gaps in the data, a method has to be developed that will compensate for the data that have no dedicated calibration measurements.

## 3. Back Calibration Method

Gould & Lyne (1998) have created average profiles for all the pulsars in the Jodrell Bank archive that have polarisation information. These integrated profiles were created by averaging all the observations of a pulsar for a certain frequency and they are assumed to be the ”true” pulse profile at that frequency for the purposes of the Back Calibration Method. Examples of such profiles can be seen in Fig. 1.

The main principle behind the Back Calibration Method is using formalisation similar to that found in Johnston (2002). For a certain time period when the telescope is observing using a certain system setup, the whole observing process introduces errors in the measurement such as those described in the previous section. To determine them, the telescope, the polarimeter and whatever electronics intervene from the point the signal enters the system to the point it is recorded are represented - to a first approximation - by a 4$`\times `$4 matrix (the so-called Mueller Matrix - M). In reality, the 16 M elements can be reduced to 8 as they are interrelated.

Therefore, the following would hold:

$$S^{}=MS$$
(1)
where $`S^{}`$ is the Stokes Parameters’ 4-vector which, when combined for all bins, will make up an observation from the same system setup period whose M is to be determined. S is the corresponding corrected (or ’true’) vector.

Thus, the M fully describing a certain system setup can be determined from Eq. 1 and a certain observation, if the ’true’ integrated profiles by Gould & Lyne (1998) are used. In particular, in order to calculate the M elements, the observation and the true profile are compared and one vector equation like Eq. 1 is formed for each bin that has a signal to noise ratio of over 7. Further equations can be formed if one compares observations from different times of day (corresponding to various parallactic angles) or from different pulsars (but always from the same system setup) with the respective true profile found in the archive.

The set of equations thus formed can then be solved for the 8 unknown elements (e.g. using a SIMPLEX method) determining in that manner M. The more equations used in the fitting process the more accurate M is likely to be. It then follows that:

$$S=M^1S^{}$$
(2)
Therefore, M can be used to retrieve the true Stokes Parameters for that particular observation from Eq. 2. The product of the inverse M times the observed Stokes Parameters will be a profile with any instrumental effects removed.

The process relies on the assumption that the archive profiles of Gould & Lyne are true representatives of the actual pulsar profiles. However, even if some of the archive profiles used are not accurate, incorporating as many as possible in the Back Calibration Method and following a selection of only the higher quality ones, ensures the required accuracy when used in conjunction with a least-squares method.

## 4. Future prospects

As soon as a procedure is fully established with which any observation existing in the database can be calibrated using information drawn from the database itself, a number of projects can then start to be examined.

For example, the precession of PSR B1828$``$11 can be verified and studied by looking at the behaviour of its PA swing through the 15-year long baseline (Stairs et al these proceedings). In addition to this, all 300 pulsars in the archive can be examined as to if they exhibit similar precession characteristics. Even outside the scope of precession, the study of the evolution of the PA swing through time (especially in such long scales) is of considerable interest on its own.

As mentioned above, the geometry of the pulsar emission can also be studied. This is feasible by fitting for the PA swings and thus determining the impact parameter $`\beta `$ (i.e. the angle between the magnetic field axis and the observer’s line of sight) and the angle $`\alpha `$ (between the magnetic axis and the rotation axis). Studying the emission geometry can provide constraints for emission models and help answer questions regarding emission regions and origin, single pulse structure and giant pulses.

Finally, it must be noted that the existing database will be extended with the use of COBRA (Coherent Online Baseband Receiver for Astronomy) the new Beowulf cluster of 182 processors. A more systematic and careful calibration is implemented from the beginning to COBRA which will allow for up to 100 MHz total B/W of coherently dedispersed and fully calibrated data with the 4 Stokes parameters of polarisation provided ’for free’ - i.e. determined in software). Such an instrument can be used in a wide number of observations where high profile resolution is necessary or to address entirely new questions such as single pulse polarimetry of Millisecond Pulsars or MSP Rotation Measures.

### Acknowledgments.

DA would like to thank Bhal Chandra Joshi and Christine Jordan for valuable discussions as well as Anita Richards without whose help the poster would not have made it to the conference. DA acknowledges the support of PPARC for the funding of this research project. IHS is supported by an NSERC UFA and Discovery Grant.

## 5. References

Gould, D. M. & Lyne, A. G. 1998 in MNRAS, 301, 1, pp. 235-260
Hamaker, J. P., Bregman, J. D., Sault, R. J. 1996 in A&AS, 117, p.137-147
Johnston, S. 2002 in Publications of the Astronomical Society of Australia, 19, 2, pp. 277-281
Stairs, I. H., Lyne, A. G., Athanasiadis D., & Kramer M. these proceedings
Stinebring, D. R., Cordes, J. M., Rankin, J. M., Weisberg, J. M., & Boriakoff, V. 1984 in A&AS, 55, June 1984, p. 247-277


text:
# Response of distance measures to the equation of state

## 1 Introduction

Current cosmological data suggests that the expansion of the universe is accelerating (e.g. Perlmutter et al. 1999; Riess et al. 1998; Efstathiou et al. 2002; Lewis & Bridle 2002; Melchiorri et al. 2002). A number of models have been proposed to explain this fact, the simplest of which is the cosmological constant. While a cosmological constant is enough to explain the current data, its constancy leads to a fine tuning problem (Sahni & Starobinski 2000; Peebles & Ratra, 2002; Padmanabhan, 2002b). A simple phenomenological generalization of the cosmological constant is to model the dark energy component that drives the acceleration as an ideal fluid with an equation of state given by $`P=w\rho `$ in which the equation of state parameter $`w`$ is allowed to vary with time. In this parameterization the cosmological constant corresponds to $`w=1`$, while for other viable models $`0.6w1`$ at the present epoch. Such a parameterization indeed arises naturally in several models, such as quintessence (Ratra & Peebles 1988; Wetterich 1988), K-essence (Armendariz-Picon et al. 1999) and a tachyonic scalar field (Gibbons 2002; Padmanabhan 2002a; Padmanabhan & Roy Choudhury 2002; Bagla et al. 2003). However the precise form of $`w(z)`$ is model dependent.

In view of this there has been considerable interest in attempts to summarize the current (and future) data in terms of a few numbers. There are several ways of doing this, such as using a polynomial approximation for $`w(z)`$ (e.g. Weller & Albrecht, 2002) or in terms of derivatives of the expansion factor (called Statefinders by Sahni et al. 2002). To be able to fit data using a polynomial form we have to choose a low order polynomial since realistic data contains only a finite amount of information. On the other hand we need to use enough polynomial orders so that the data is adequately described. Saini et al. (2003) investigate in some detail the issue of how to go about this and show that the current data and near future data seem to require at most a low order polynomial. (This, of course, does not imply that the true equation of state is also of a low order polynomial form.)

These methods attempt to capture the effect of an unknown function $`w(z)`$, which is equivalent to infinite number of parameters, in terms of a finite (and often small) number of parameters. In the limited redshift range probed by supernova data the true variations in $`w(z)`$ might be such that a low order polynomial fit is a reasonable approximation. Alternatively one can expand the function $`w(z)`$ in terms of a some set of basis functions which are complete in the given redshift interval. If the basis functions are chosen judiciously only a small set of expansion coefficients will be required to describe the function $`w(z)`$ within the limits of the experimental accuracy. Owing to the random noise in cosmological data any of these choices could fit the data adequately in a maximum likelihood sense. Since the distance measures are integrals over nonlinear expressions involving $`w(z)`$, any parameter determined by such a maximum likelihood analysis represents some non-trivial average of the true $`w(z)`$. Our main aim in this paper is to show that the actual *functional* relation of the cosmological distance measures to $`w(z)`$ is not far from being linear. This allows us to deduce an approximate relation between the fitted coefficients and the true $`w(z)`$. This allows one to obtain a quantitative description of the averaging involved in reducing a function $`w(z)`$ to a finite set of numbers.

The plan of this paper is as follows. In Section 2 we derive the linear response approximation relating the coordinate distance and the equation of state of the dark energy, and show that within a reasonable range of parameters for the dark energy it works to better than $`2\%`$. In Section 3 we use this approximation to relate the fitted coefficients of a polynomial form to the true equation of state for the case of fitting to the luminosity distance, and show that they are related through a weighted integral of the true equation of state. We then generalize this result to the case of an arbitrary functional form which is linear in the parameters. In Section 4 we derive a similar relationship for CMB where, in the simplest case, only one measurement of the angular distance to the last scattering surface is available. In Section 5 we show that the linear approximation also enables a non-parametric estimation of $`w(z)`$. Our conclusions are presented in Section 6.

## 2 Response of geometry to the equation of state

The luminosity and angular diameter distances depend on the dark energy through the coordinate distance $`r[w,z]`$ according to the equations: $`D_L(z)=(1+z)r[w,z]`$ and $`D_A(z)=r[w,z]/(1+z)`$. The square brackets in this notation explicitly shows that at any redshift $`z`$ the coordinate distance requires full knowledge of the function $`w(z)`$ In a flat universe the coordinate distance $`r[w,z]`$ is given by

$$r[w,z]=(1+g)^{1/2}_1^{1+z}𝑑xx^{3/2}[g+Q[w,x]]^{1/2}$$
(1)
where $`x=1+z`$, $`g=\mathrm{\Omega }_\mathrm{m}/(1\mathrm{\Omega }_\mathrm{m})`$,

$$Q[w,z]=\mathrm{exp}[3_1^{1+z}𝑑xw(x)/x],$$
(2)
and we have set $`c`$ and $`H_0`$ equal to unity. To explore the behaviour of $`r[w,z]`$ when different equation of state functions $`w(z)`$ are used, we need to understand the sensitivity of $`r`$ to $`w(z)`$. This can be characterized by the *functional derivative* of $`r`$ with respect to $`w(z)`$. Since this is not a routine weapon in the arsenal of the astronomer, we shall briefly introduce the concept before applying it.

In the case of a real function $`f(x)`$, the sensitivity of the function to the independent variable $`x`$ can be characterized by the derivative $`df/dx=f^{}(x)`$. Broadly speaking, a large value for $`f^{}(a)`$ indicates that $`f`$ is relatively more sensitive to the independent variable around the point $`x=a`$; and a small value for $`f^{}(a)`$ will indicate relative insensitivity of $`f`$ to the independent variable around $`x=a`$. This follows directly from the definition of derivative of a function

$$f^{}(a)=\underset{ϵ0}{lim}\frac{f(a+ϵ)f(a)}{ϵ}.$$
(3)
In the case of $`f`$ depending not on a single variable but on a *function* $`p(x)`$ we need to study how $`f`$ changes if the function $`p(x)`$ is changed slightly around a point $`x=b`$. This is best done by changing the function $`p(x)`$ to a new function $`p_1(x)p(x)+ϵ\delta _\mathrm{D}(xb)`$ which adds a “spike” at $`x=b`$ with a strength proportional to $`ϵ`$. We can now evaluate the value of $`f`$ for both $`p(x)`$ and $`p_1(x)`$. The difference in the numerical values of $`f`$ in the limit of $`ϵ0`$ is a good measure of the sensitivity of $`f`$ to the functional form of $`p(x)`$ around $`x=b`$. More formally, this functional derivative is defined by

$$\frac{\delta f[p,x]}{\delta p(b)}=\underset{ϵ0}{lim}\frac{f[p+ϵ\delta _\mathrm{D}(xb),x]f[p,x]}{ϵ}.$$
(4)
Just as the ordinary derivative of a function depends on the location at which it is evaluated, the functional derivative depends both on the form of $`p(x)`$ around which the sensitivity is measured, $`x=a`$, as well as the point $`x=b`$ at which the input function is perturbed.

For the purposes of this paper we need the response of the coordinate distance $`r[w,z]`$ at a redshift $`z`$ to a change in the equation of state at a different redshift, $`z^{}`$. This can be computed around a given fiducial $`w(z)=w^{\mathrm{fid}}(z)`$ from the functional derivative

$$\frac{\delta r[w^{\mathrm{fid}},z]}{\delta w(z^{})}=\underset{ϵ0}{lim}\frac{r[w^{\mathrm{fid}}+ϵ\delta _\mathrm{D}(zz^{}),z]r[w^{\mathrm{fid}},z]}{ϵ}$$
(5)
defined exactly as in equation (4) (also see Huterer & Turner, 2001 for a similar application of this idea)

For the rest of this section we switch from discussing the coordinate distance to the luminosity distance and change the independent variable from redshift $`z`$ to $`x=1+z`$. Multiplying by $`(1+z)`$ we obtain the response function for the luminosity distance $`\delta D_L[w^{\mathrm{fid}},x]/\delta w(x^{})K_w(x,x^{})`$. The subscript $`w`$ on the kernel denotes that it depends on $`w^{\mathrm{fid}}`$, the fiducial equation of state around which the approximation holds. Evaluating the expression in Eq 5 we obtain

$$K_w(x,x^{})=\{\begin{array}{c}\frac{3x(1+g)^{1/2}}{2x^{}}_x^{}^x\frac{dy}{y^{3/2}}\frac{Q[w^{\mathrm{fid}},y]}{\left(g+Q[w^{\mathrm{fid}},y]\right)^{3/2}}\hfill \\ 0(\text{for}x<x^{})\hfill \end{array}.$$
(6)
The kernel $`K_w`$ is a function of two parameters, the redshift at which $`w(z)`$ is perturbed and the redshift at which we consider the change in the luminosity distance. A surface plot of this function is shown in Fig 1. Since the effect of varying $`w(z)`$ at a redshift $`z^{}`$ is felt only at $`z>z^{}`$,the kernel is zero in half the plane. For small $`\delta w`$ we can use this result to approximate the calculation of the luminosity distance as

$`D_L[w^{\mathrm{fid}}+\delta w,z]`$ $``$ $`D_L[w^{\mathrm{fid}},z]+\delta D_L`$
$``$ $`D_L[w^{\mathrm{fid}},z]+{\displaystyle _0^z}K_w(z,z^{})\delta w(z^{})𝑑z^{}.`$
To be useful this linear response approximation should hold to a good accuracy for a reasonable range of $`w(z)`$. The SuperNova Acceleration Probe (SNAP) survey is expected to observe about $`2000`$ Type 1a supernovae (SNe), up to a redshift $`z1.7`$, each year (Aldering et al. 2002). A single supernova will measure the luminosity distance with a relative error of $`7\%`$. If we bin the supernovae in redshifts interval of $`0.02`$, this will give a relative error in the luminosity distance of about $`1\%`$. Saini et al. (2002) show that given this level of precision and given the present uncertainty in the value of $`\mathrm{\Omega }_m`$, the data seems to require at the most a linear polynomial order in $`w(z)`$. To show how accurately Eq 2 gives the luminosity distance in this restricted case, in Fig 2 we plot in the $`w_0`$-$`w_1`$ plane the percentage difference between the true luminosity distance and the one computed through Eq 2 at $`z=1.5`$. The kernel $`K`$ was computed for $`w^{\mathrm{fid}}(z)=0.5`$ in this calculation. Within the region $`1<w(z)<0`$ the approximation is better than 2%. Therefore, the accuracy of the linear approximation is quite acceptable with respect to the projected uncertainties in the luminosity distance obtained from a SNAP class experiment. The approximation works even better at smaller redshifts but less so at higher redshifts, but the measurement errors are expected to be larger here. Therefore, the conclusions drawn from this approximation are expected to remain valid until an unprecedented accuracy is achieved in the measurement of luminosity distance.

## 3 Interpretation of the fitting parameters

Since we have no fundamental understanding of the nature of dark energy, it is necessary to assume some suitable, versatile, functional form for $`w(z)`$ to fit the cosmological data. Since the luminosity distance depends on the equation of state through an integral relation, and the cosmological data is often fitted through maximum likelihood method, we can always add a small fast varying term to any well fitting $`w(x)`$, while still retaining a good fit 111 For example $`w(x)`$ and $`w(x)+\theta \mathrm{sin}(x/\theta )`$, where $`\theta 0`$, would both fit the data equally well though none of the derivatives agree. Although this is a contrived example this point has been well illustrated more realistically in Maor et al. (2000).. Due to this integral dependence any flexible parameterization eventually recovers only some integrated property of the underlying equation of state. The simple polynomial forms for $`w(z)`$ usually assumed for the purposes of fitting the cosmological data are often viewed as the first few terms in a series expansion of $`w(z)`$. Since the behaviour of $`w(z)`$ is not necessarily polynomial like, the recovered coefficients of the polynomial are not related in a simple manner to the true coefficients of the series expansion. In this Section we use the linear response approximation described above to compute the expectation value for the coefficients of the polynomial approximating the $`w(z)`$, and relate them to the “true” input $`w(z)`$. This will show that such fitting functions serve a useful purpose of measuring some integrated properties of the true equation of state of the dark energy.

### 3.1 Interpreting the fit with a constant $`w`$

The simplest fit to the cosmological data is the *constant* $`w`$ model, $`w^{\mathrm{fit}}(z)=w_0^{\mathrm{fit}}`$. We are interested in relating this to the true $`w(z)=w^{\mathrm{true}}(z)`$. For an arbitrary equation of state this relation is non-trivial and non-linear but given the approximations considered in § 2, we can construct the relationship analytically. Let us first consider fitting the luminosity distance to redshift. This requires a knowledge of the present day matter density $`\mathrm{\Omega }_m`$ as well and, for simplicity, we shall assume that this is known. (When $`\mathrm{\Omega }_m`$ is unknown, both $`\mathrm{\Omega }_m`$ and $`w`$ become biased; see Maor et al. 2001). Our approach could be extended to the case in which $`\mathrm{\Omega }_\mathrm{m}`$ is not known, however this is beyond the scope of the present paper.

In the maximum likelihood reconstruction the quantity

$`\chi ^2={\displaystyle _0^{z_{\mathrm{max}}}}𝑑z\left({\displaystyle \frac{D_L^{\mathrm{fit}}(z)(D_L^{\mathrm{true}}(z)+n(z))}{\sigma (z)}}\right)^2`$ (8)
is minimized, where $`n(z)`$ is the noise in the measurement at redshift $`z`$ and $`\sigma (z)`$ is the variance. (We have replaced the conventional summation with an integral.) We now approximate the fitting function and the luminosity distance by

$`D_L^{\mathrm{fit}}(z)`$ $``$ $`D_L[w^{\mathrm{fid}},z]+\mathrm{\Delta }w{\displaystyle _0^{z_{\mathrm{max}}}}K_w(z,z^{})𝑑z^{}`$ (9)
$`D_L(z)`$ $``$ $`D_L[w^{\mathrm{fid}},z]+{\displaystyle _0^{z_{\mathrm{max}}}}K_w(z,z^{\prime \prime })\delta w(z^{\prime \prime })𝑑z^{\prime \prime }`$ (10)
where we have assumed a constant $`w^{\mathrm{fid}}`$ and defined $`\mathrm{\Delta }w=w_0^{\mathrm{fit}}w^{\mathrm{fid}}`$ and $`\delta w(z)=w^{\mathrm{true}}(z)w^{\mathrm{fid}}`$. Minimizing $`\chi ^2`$ and taking the expectation value we obtain

$`\mathrm{\Delta }w={\displaystyle _0^{z_{\mathrm{max}}}}\mathrm{\Phi }_w(z^{})\delta w(z^{})𝑑z^{}`$ (11)
where

$`\mathrm{\Phi }_w(z)={\displaystyle \frac{K_w(x,x^{})K_w(x,1+z)/\sigma ^2(x)𝑑x𝑑x^{}}{K_w(x,x^{})K_w(x,x^{\prime \prime })/\sigma ^2(x)𝑑x𝑑x^{}𝑑x^{\prime \prime }}}`$ (12)
where all integrals are in the range $`(0,z_{\mathrm{max}})`$. The noise term has dropped out since its expectation value is zero. Adding $`w_0^{\mathrm{fit}}`$ to both sides of Eq. (11) we obtain

$`w_0^{\mathrm{fit}}={\displaystyle _0^{z_{\mathrm{max}}}}\mathrm{\Phi }_w(z^{})w^{\mathrm{true}}(z^{})𝑑z^{}.`$ (13)
The fitted constant, $`w^{\mathrm{fit}}`$, is just a weighted average of the true equation of state, with weighting function $`\mathrm{\Phi }_w`$. This weighting function is shown in Fig. 3 for the case of supernovae distributed evenly from $`z=0`$ to $`z=2.0`$. It decreases steadily with redshift, as might be expected from the fact that all the supernovae are affected by the equation of state at $`z=0`$, whereas the equation of state at $`z>2.0`$ affects no supernovae. This result shows that the value of $`w_0^{\mathrm{fit}}`$ obtained by fitting to supernovae data remains invariant if the $`w^{\mathrm{true}}(z^{})`$ is changed in such a way that the weighted integral is unaffected.

If the linear response approximation did not hold to a high accuracy then Eq. (13) would lead to different results depending on the value for $`w^{\mathrm{fid}}`$ used to compute the weighting function $`\mathrm{\Phi }_w`$. Therefore a further test for the linear response approximation is provided by the variation of the weighting function evaluated for a different $`w^{\mathrm{fid}}`$. In Fig 3 we show the weighting function computed for two different values of $`w^{\mathrm{fid}}`$, namely $`w^{\mathrm{fid}}=1.0`$ and $`w^{\mathrm{fid}}=0`$. The figure shows that the variation in the weighting function $`\mathrm{\Phi }_w`$ is small.

To further quantify this difference we consider the change in $`w_0^{\mathrm{fit}}`$ with the change in the weighting function

$`\mathrm{\Delta }w_0^{\mathrm{fit}}={\displaystyle _0^{z_{\mathrm{max}}}}\mathrm{\Delta }\mathrm{\Phi }_w(z^{})w^{\mathrm{true}}(z^{})𝑑z^{}.`$ (14)
Squaring both sides and using the Schwartz inequality we obtain the bound:

$`(\mathrm{\Delta }w_0^{\mathrm{fit}})^2{\displaystyle _0^{z_{\mathrm{max}}}}(\mathrm{\Delta }\mathrm{\Phi }_w)^2𝑑z{\displaystyle _0^{z_{\mathrm{max}}}}(w^{\mathrm{true}})^2𝑑z`$ (15)
Since we have assumed that $`|w|<1`$ in the relevant range, the second integral on the right hand side has a maximum value $`z_{\mathrm{max}}`$. (We have computed the weighting function for $`z_{\mathrm{max}}=2.0`$.) With these numbers we can now compute $`(\mathrm{\Delta }w^{\mathrm{fit}})^2`$ by considering the difference of $`\mathrm{\Phi }_w`$ evaluated at $`w_0=1`$ and $`w_0=0`$. We obtain $`|\mathrm{\Delta }w^{\mathrm{fit}}|0.1`$. This shows that the expected error in the determination of $`w`$ through the linear response approximation is about $`10\%`$ if we linearize the luminosity distance function around an *incorrect* initial guess. In practice, of course, one could do much better by first fitting the cosmological data with a constant $`w`$ and then expanding around that value so that the approximation works better than when it is expanded about an arbitrary point.

### 3.2 Generalization

We have derived above an expression for the expectation value of $`w^{\mathrm{fit}}`$ for the case of a *constant* $`w`$ fit. The corresponding expressions for the case of more complicated fitting functions are, however, quite cumbersome. We can simplify the notation a little by considering discretized expressions. Suppose the luminosity distance is known at a large number of, uniformly distributed redshifts $`z_i`$ where $`i=1,N`$. We consider a model $`w(z)`$ which is given at redshifts $`z_k^{}`$ where $`k=1,M`$.

After having chosen a $`w^{\mathrm{fid}}`$, which for simplicity is taken to be a constant, we can define the vectors $`𝐝\{\delta D_L(z_i)\}`$ and $`𝐰\{\delta w(z_i^{})\}`$ and matrix $`𝒦\{K(z_i,z_j^{})\delta z\}`$, where $`𝒦`$ is a $`N\times M`$ matrix and $`\delta z`$ denotes the redshift interval between the bins. With these definitions Eq (2) becomes

$$𝐝=𝒦𝐰,$$
(16)
where both $`𝐝`$ and $`𝐰`$ are small quantities. Next we give the equivalent discrete version of the weighting function obtained in the previous subsection. We have $`w^{\mathrm{fit}}(z)=w^{\mathrm{fid}}+\mathrm{\Delta }w`$, where $`\mathrm{\Delta }w`$ is a constant; therefore, using the notation defined above we have $`𝐰^{\mathrm{fit}}=𝐮\mathrm{\Delta }w`$, where $`𝐮=(1,1,\mathrm{},1)`$, and similarly $`w^{\mathrm{true}}=w^{\mathrm{fid}}+\delta w(z)`$, $`𝐰^{\mathrm{true}}=\{\delta w(z_i)\}`$. With these definitions Eq. (11) and (12) translate to

$`\mathrm{\Delta }w={\displaystyle \frac{𝐮^\mathrm{T}}{𝐮^\mathrm{T}𝐮}}𝐰^{\mathrm{true}};=𝒦_w^\mathrm{T}𝒦_w.`$ (17)
If instead we fit the data with a linear model $`w^{\mathrm{fit}}=w^{\mathrm{fid}}+\mathrm{\Delta }w_0+\mathrm{\Delta }w_1z`$ then, following the same procedure that led to Eq. 11, we find

$`\mathrm{\Delta }w_0`$ $`=`$ $`{\displaystyle \frac{\left[(𝐳^\mathrm{T}𝐳)𝐮^\mathrm{T}(𝐮^\mathrm{T}𝐳)𝐳^\mathrm{T}\right]}{(𝐮^\mathrm{T}𝐮)(𝐳^\mathrm{T}𝐳)(𝐮^\mathrm{T}𝐳)^2}}𝐰^{\mathrm{true}}`$ (18)
$`\mathrm{\Delta }w_1`$ $`=`$ $`{\displaystyle \frac{\left[(𝐮^\mathrm{T}𝐮)𝐳^\mathrm{T}(𝐮^\mathrm{T}𝐳)𝐮^\mathrm{T}\right]}{(𝐮^\mathrm{T}𝐮)(𝐳^\mathrm{T}𝐳)(𝐮^\mathrm{T}𝐳)^2}}𝐰^{\mathrm{true}},`$ (19)
where we have also defined the vector $`𝐳=\{z_i\}`$. This generalizes the concept of the weighting function to the linear case. It should be noted that the weighting function for $`\mathrm{\Delta }w_0`$ is different than in Eq. (17), since we are also fitting for $`w_1`$. However, as expected, these results agree for the special case when $`𝐰^{\mathrm{true}}=\alpha 𝐮`$.

In the most general case we could approximate $`w(z)`$ as a linear combination of an arbitrary set of functions $`F_i`$ (e.g. Gerke & Efstathiou, 2002) as

$`w^{\mathrm{fit}}(z)=w^{\mathrm{fid}}+{\displaystyle \underset{i=1}{\overset{N}{}}}c_iF_i(z)`$ (20)
In this case the coefficients $`c_i`$s are related to the true equation of state as follows. Defining the $`N`$ vectors $`𝐅_i=\{F_i(z_k)\}`$ we obtain $`𝐰^{\mathrm{fit}}=_{i=1}^Nc_i𝐅_i`$. On performing a maximum likelihood analysis we find an expression for the coefficients $`c_i`$ as follows

$`c_i={\displaystyle \underset{i=1}{\overset{N}{}}}a_{ij}^1y_j`$ (21)
where we have defined the two matrices

$`a_{ij}`$ $`=`$ $`𝐅_i^\mathrm{T}𝐅_j`$ (22)
$`y_i`$ $`=`$ $`𝐅_i^\mathrm{T}𝐰^{\mathrm{true}}.`$ (23)
## 4 Effective $`w`$ seen by the CMB

Current and near future supernova surveys will probe redshifts only up to around $`z=2`$. We can also compute the effective equation of state probed by the Cosmic Microwave Background (CMB) at $`z_{\mathrm{cmb}}1000`$. In the simplest case the CMB gives us the angular distance to the last scattering surface and therefore it does not give detailed information on the behaviour of $`w(z)`$. Huey et al. (1999) show that a $`w_{\mathrm{eff}}`$ that keeps the CMB temperature and the matter power spectrum unaffected to an accuracy of $`5`$$`10\%`$, is given by

$$w_{\mathrm{eff}}=\frac{\mathrm{\Omega }_Q(a)w(a)𝑑a}{\mathrm{\Omega }_Q(a)𝑑a}$$
(24)
where the integral extends up to the CMB redshift. This is in the same spirit as the weighting function corresponding to fitting a constant constant equation of state described earlier. Following a similar procedure to that in the previous section we obtain

$`w_{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \mathrm{\Phi }^{\mathrm{cmb}}(z)w(z)𝑑z}`$ (25)
$`\mathrm{\Phi }^{\mathrm{cmb}}(z)`$ $`=`$ $`{\displaystyle \frac{K^{}(1+z_{\mathrm{cmb}},1+z)}{_0^{z_{\mathrm{cmb}}}K^{}(1+z_{\mathrm{cmb}},1+z)𝑑z}}`$ (26)
where $`K^{}=K/(1+z)^2`$ due to the difference between $`D_L`$ and $`D_A`$.

In Figure 4 (solid line) we plot this weighting function as calculated using the linear approximation by expanding around $`w=1`$. It has the same rough shape as the weighting function for supernovae but this time extends to higher redshift. We see that despite the high redshift of the CMB the weighting function falls off close to zero by redshift $`z5`$. This is due to the fact that the dark energy term in the Friedman equations is increasingly unimportant as we go to high redshifts. For example it is well known that the cosmological constant is dynamically unimportant at CMB redshifts.

Overlaid on Fig. 4 (dashed line) is the weighting function from Eq. 24, converting the integral from the scale factor $`a`$ to redshift $`z`$. The good agreement shows that weighting with $`\mathrm{\Omega }_Q`$ works well for models close to the cosmological constant. In Fig 5 we show the same plot but this time expanding around $`w^{\mathrm{fid}}=0`$. The two results do not agree to a good accuracy. This shows that the weighting with $`\mathrm{\Omega }_Q`$ should not work so well for models which have an average $`w`$ closer to zero. These figures also show that the linear response approximation works less well for the CMB since the weighting function changes by a large amount when we expand around $`w=1`$ and $`w=0`$. This means that the weighting function would only apply for those models for which $`|ww_{\mathrm{eff}}|`$ is small.

## 5 Non-parametric reconstruction of the equation of state parameter

Although the fitting functions that we have considered so far provide useful information about the true equation of state of the dark energy, the linear approximation discussed above in principle allows a non-parametric reconstruction of the equation of state.

Consider a given SNe data set. To a first approximation we can fit it with a constant $`w`$, even though this may not be a good fit. We set the fiducial equation of state $`w^{\mathrm{fid}}`$ to this best fit constant value in all that follows. We next compute the difference between the SNe distances and the $`D_L`$ computed for $`w^{\mathrm{fid}}`$ to obtain the residuals $`\delta D_L`$. We use the notation described in the first paragraph of § 3.2 so these are given by $`𝐝^s`$, where the superscript $`s`$ signifies that these are obtained from measured supernovae distances. These are related to the $`w(z)`$ in bins through Eq. 16, where $`𝐰=\{w(z_k)w_0\}`$ with $`w_0`$ given by the constant $`w`$ fit. Since these are noisy estimates of distances we need to minimize the $`\chi ^2`$ function with respect to $`𝐰`$ to obtain the Maximum Likelihood estimator for $`𝐰`$; that is, we minimize $`\chi ^2=(𝐝^s𝒦𝐰)^T(𝐝^s𝒦𝐰)`$ and obtain

$`𝐰=(𝒦^\mathrm{T}𝒦)^1𝒦^\mathrm{T}𝐝^s,`$ (27)
where we have assumed the noise on $`𝐝^s`$ is constant with redshift, although the result could be extended for the general case.

This equation also allows the trivial calculation of the Fisher matrix $`=𝒦^\mathrm{T}𝒦`$ corresponding to the uncertainties on the reconstructed equation of state. It was noted by Huterer & Starkman (2002) that the Fisher matrix was a surprisingly weak function of the model parameters which matches with the result obtained above and shows that it reflects on the validity of the linear response approximation.

Although Eq 27 gives a formal solution of the problem, this estimation is, in general, very noisy. Huterer et al. show that even for SNAP-like data there are only a few principal components which are well determined by this method. It is clear that this lack of resolution is largely due to the fact that no constraints are imposed on the behaviour of $`𝐰`$. We discuss some of the ways of rectifying this in a separate paper.

## 6 Conclusions

We have shown that in the relevant range of parameters expected for the dark energy the luminosity distance is a linear functional to the equation of state to a surprisingly high level of accuracy. This approximation allows us to find the relationship between the usual polynomial models for $`w(z)`$ and the true underlying equation of state of the dark energy. Although the usual interpretation of these polynomial approximation is in terms of a series expansion of $`w(z)`$, we show that the coefficients of the polynomial approximation are related to the true equation of state through an integral relation. Only in the exact polynomial-like equation of state would these approximations measure the real $`w(z)`$. We show that the fitting of cosmological data with such forms is still useful since the parameters of such models measure certain, well-defined, integrated properties of the underlying equation of state. Finally, the approximation allows a formal, non-parametric way to measure the equation of state.

### ACKNOWLEDGMENTS

TDS and SLB acknowledges financial support from PPARC. SLB also thanks Selwyn College for support in the form of a Trevelyan Research Fellowship. We thank Irit Maor, Jochen Weller and Shiv Sethi for helpful discussions.


text:
# Entanglement Entropy and The Density Matrix Renormalization Group

## 1 Entanglement entropy in some relevant geometries

Entanglement or nonseparability refers to the existence of quantum correlations between two sets of degrees of freedom of a physical system that can be considered as subsystems. It is natural that two (sub)systems in interaction are entangled an that they still are entangled after their interaction has ceased. Particularly interesting situations arise when two entangled systems become causally disconnected because of an event horizon.

### 1.1 Introduction: entanglement entropy

The entanglement of two parts of a quantum system can be measured by the von Neumann entropy. This is defined in terms of the density matrices of either part. Let us consider, for later convenience, one part as “left” or “interior” and another as “right” or “exterior”, in a yet imprecise sense. Then, let us represent states belonging to the left (or interior) with small letters and states belonging to the right (or exterior) with capital letters. A basis for the global states (left plus right) is $`\{|a\}\{|A\}`$. By representing the ground state in this basis as

$$|0=\underset{aA}{}\psi _{aA}|a|A,$$
(1)
we define a coefficient matrix $`\psi _{aA}`$. Then we have two different density matrices:

$$\rho _R=\frac{\psi ^{}\psi }{\mathrm{Tr}\psi ^{}\psi },\rho _L=\frac{\psi ^{}\psi ^^T}{\mathrm{Tr}\psi ^{}\psi ^^T}.$$
(2)
Correspondingly, we have two von Neumann entropies:

1. $`S_R=\mathrm{Tr}_A\left(\rho _R\mathrm{ln}\rho _R\right)`$
2. $`S_L=\mathrm{Tr}_a\left(\rho _L\mathrm{ln}\rho _L\right)`$

Now it is important to recall the “symmetry theorem”, which states that both entropies are equal, $`S_R=S_L`$ (this can be proved in several ways, the most popular one appealing to the Schmidt decomposition of the entangled state. The equality of entropies implies that they are associated with properties shared by both parts, that is, with (quantum) correlations.

More generally, two non-interacting parts of a quantum system can be originally in respective mixed states. After their interaction, which we describe as an arbitrary unitary evolution of the composite system, the initial density matrix $`\rho _L\rho _R`$ has evolved to $`\rho _{LR}^{}`$. It is easy to see that the partial traces $`\rho _L^{}`$ and $`\rho _R^{}`$, in general, have von Neumann entropies $`S_L^{}`$ and $`S_R^{}`$ such that $`S_L^{}+S_R^{}S_L+S_R`$ . Of course, if the initial state is pure $`S_L=S_R=0.`$

### 1.2 Field theory half-space density matrix

Let us now consider the quantum system to be a chain of coupled oscillators. Moreover, we shall chiefly work in the continuum limit, where the concepts and mathematical expressions are more transparent, in spite of dealing with non-denumerable sets of degrees of freedom (we will return to a discrete chain to describe the density matrix renormalization group algorithm). The action for this model, namely, a one-dimensional scalar field, is

$$A[\phi (x,t)]=𝑑t𝑑x\left(\frac{1}{2}\left[(_t\phi )^2(_x\phi )^2\right]V(\phi )\right),$$
(3)
where $`\phi `$ is the field.

Let us obtain a path integral representation for the density matrix on the half-line of a system that is in its ground state . In the continuum limit, the half-line density matrix is a functional integral,

$$\rho [\phi _R(x),\phi _R^{}(x)]=D\phi _L(x)\psi _0[\phi _L(x),\phi _R(x)]\psi _0^{}[\phi _L(x),\phi _R^{}(x)],$$
(4)
where the subscripts refer to the left or right position of the coordinates with respect to the boundary (the origin). Now, we must express the ground-state wave-function as a path integral,

$$\psi _0[\phi _L(x),\phi _R(x)]=D\phi (x,t)\mathrm{exp}\left(A[\phi (x,t)]\right),$$
(5)
where $`t(\mathrm{},0]`$ and with boundary conditions $`\phi (x,0)=\phi _L(x)`$ if $`x<0`$, and $`\phi (x,0)=\phi _R(x)`$ if $`x>0`$. The conjugate wave function is given by the same path integral and boundary conditions but with $`t[0,\mathrm{})`$. Substituting into Eq. (4) and performing the integral over $`\phi _L(x)`$, one can express $`\rho (\phi _R,\phi _R^{})`$ as a path integral over $`\phi (x,t)`$, with $`t(\mathrm{},\mathrm{})`$, and boundary conditions $`\phi _R(x,0+)=\phi _R^{}(x)`$, $`\phi _R(x,0)=\phi _R(x)`$. In other words, $`\rho (\phi _R,\phi _R^{})`$ is represented by a single path integral covering the entire plane with a cut along the positive semiaxis, where the boundary conditions are imposed.

### 1.3 Angular quantization and Rindler space

In Euclidean two-dimensional field theory, the generator of rotations in the $`(x,t)`$ plane is given by

$$=𝑑x(tT_{11}xT_{00}),$$
(6)
in terms of the components of the stress tensor computed from the action (3). To simplify, one can evaluate it at $`t=0`$. Let us consider a free action ($`V=0`$). In the Schrödinger representation, we should replace the momentum $`\mathrm{\Pi }=_t\phi `$ with $`\mathrm{\Pi }(x)=i\delta /\delta \phi (x)`$. However, as in canonical quantization, one rather uses the second-quantization method, which diagonalizes the Hamiltonian by solving the classical equations of motion and quantizing the corresponding normal modes. Let us recall that, in canonical quantization, if we disregard anharmonic terms, the classical equations of motion in the continuum limit become the Klein-Gordon field equation, giving rise to the usual Fock space. Not surprisingly, the eigenvalue equation for $``$ leads to the Klein-Gordon equation in polar coordinates in the $`(x,t)`$ plane,

The free field wave equation in polar coordinates,

$$(\mathrm{\Delta }+m^2)\phi =\left(\frac{1}{r}\frac{}{r}r\frac{}{r}+\frac{1}{r^2}\frac{^2}{\varphi ^2}+m^2\right)\phi =0,$$
(7)
can be solved by separating the angular variable: it becomes a Bessel differential equation in the $`r`$ coordinate with complex solutions $`I_{\pm i\mathrm{}}(mr)`$, $`\mathrm{}`$ being the angular frequency. We have a continuous spectrum, which becomes discrete on introducing boundary conditions. One of them must be set at a short distance from the origin, to act as an ultraviolet regulator , necessary in the continuum limit.

Therefore, the second-quantized field is (on the positive semiaxis $`t=0\varphi =0`$, $`xr`$)

$$\phi (x)=\frac{d\mathrm{}}{2\pi }\frac{b_{\mathrm{}}I_i\mathrm{}(mx)+b_{\mathrm{}}^{}I_i\mathrm{}(mx)}{\sqrt{2\mathrm{sinh}(\pi \mathrm{})}},$$
(8)
where we have introduced annihilation and creation operators and where the term that appears in the denominator is just for normalization, to ensure that those operators satisfy canonical conmutations relations. There is an associated Fock space built by acting with $`b_{\mathrm{}}^{}`$ on the “vacuum state”. These states constitute the spectrum of eigenstates of $``$, which adopts the form $`=𝑑\mathrm{}\mathrm{}b_{\mathrm{}}^{}b_{\mathrm{}}`$ (where the integral is replaced with a sum for discrete $`\mathrm{}`$).

The type of quantization just exposed was first introduced in the context of quantization in curved space, in particular, in Rindler space. Rindler space is just Minkowski space (therefore, not curved) in coordinates such that the time is the proper time of a set of accelerated observers. Its interesting feature is the appearance of an event horizon, which implies that the ground state (the Minkowski vacuum) is a mixed (thermal) state . The connection with black hole entropy and Hawking radiation is explained in the next section.

It is pertinent to note that the functions $`I_{\pm i\mathrm{}}(mx)`$ have wave-lengths that increase with $`x`$. It is illustrative to represent a real “angular wave”,

$$K_i\mathrm{}(mx)=\frac{i\pi }{2\mathrm{sinh}(\pi \mathrm{})}[I_i\mathrm{}(mx)I_i\mathrm{}(mx)].$$

This solution is oscillatory for $`x<\mathrm{}/m`$, with a wavelength proportional to $`x`$, and decays exponentially for $`x>\mathrm{}/m`$ (Fig. 1). The fact that the wave length vanishes at $`x=0`$ is to be expected from the Rindler space viewpoint, because it corresponds to the horizon.

### 1.4 Black hole entropy

Let us consider the Schwarzschild geometry in the Kruskal-Szekeres coordinates $`u,v`$, defined by

$$uv=16M^2\left(\frac{r}{2M}1\right)\mathrm{exp}\left(\frac{r}{2M}1\right),\frac{u}{v}=\mathrm{exp}\frac{t}{2M}.$$

The horizon is given by $`u=0`$ or $`v=0`$, like in the Rindler geometry of the previous section. Then, if we further define

$$Z=u+v,T=uv,$$

these coordinates behave like the ordinary Minkowskian coordinates. Moreover, for small $`u`$ or $`v`$ (or the large $`M`$ limit), the curvature can be neglected and the geometry becomes locally the one of Rindler space. Therefore, we could define radial coordinates $`Z=\rho \mathrm{cosh}\tau ,`$ $`T=\rho \mathrm{sinh}\tau ,`$ and perform a radial quantization like we did in the preceding section.

Once established that the geometry near the black-hole horizon is locally the Rindler geometry of the preceding section, we can readily transfer the form of the density matrix of a scalar field therein, where we now ignore (trace over) the degrees of freedom inside the horizon. Hence, we can define a von Neumann entropy associated with this density matrix. Furthermore, in so doing, we can appreciate that the concept of black-hole entropy takes a new meaning: in addition to being of quantum origin, this entropy is related with shared properties between interior and exterior, namely, with the horizon. In addition, the radial vacuum is a thermal state with respect to the original Schwarzshild coordinates, giving rise to Hawking radiation .

## 2 Quantum information and RG transformations

### 2.1 Information theory and maximum entropy principle

The entropy concept appeared in Thermodynamics but only took a truly fundamental meaning with the advent of information theory. In this theory, entropy is just missing information, while information itself is often called negentropy. To recall basic definitions, the information attached to an event that occurs with probability $`p_n`$ is $`I_n=\mathrm{ln}p_n.`$ Hence, the average information (per event) of a source of events is $`S(\{p_n\})=_np_nI_n=_np_n\mathrm{ln}p_n.`$

The previous definitions, given by Shannon in his theory of communication, seem unrelated to thermodynamic entropy as a property of a physical system. However, according to the foundations of Statistical Mechanics on Probability Theory (the Gibbs concept of ensembles), a clear relation can be established, as done by Jaynes . Jaynes made connection with the Bayesian philosophy of probability theory, in which the concept of “a priory” knowledge is crucial. Indeed, although the exact microscopic state of a system with many degrees of freedom may be unknown, one has some “a priory” knowledge given by the known macroscopic variables. This is a particular case of Jaynes’ adaptation of the Bayesian probability theory, wich postulates that the best probability distribution to be attributed to a stochastic event is such that it incorporates only the “a priory” knowledge about the event and nothing else. This postulate amounts to Jaynes’ maximum entropy principle: given some constraints, one must find the maximum entropy probability distribution (density matrix, in the quantum case) compatible with those constraints, usually, by implementing them via Lagrange multipliers. In particular, more constraints mean less missing information and lead to less entropy.

### 2.2 Quantum information

The concepts of Shannon’s classical theory of communication have quantum analogues . Nevertheless, the quantum theory of communication is richer (and less intuitive!). Indeed, the key new notion in the quantum theory is entanglement (already described in the foregoing): if a state (an event) is entangled with the environment, we have the type of purely quantum phenomena to which the EPR paradox is associated.

Schumacher posed the problem of communication of an entangled state . (The technical name is transposition, since the copy of a quantum state is not possible: no-cloning theorem .) His conclusion is that the von Neumann entropy of the state is the quantity that determines the fidelity of the transposition: it is possible to transpose the state with near-perfect fidelity if the signal can carry at least that information.

The fidelity is simply defined as the overlap between normalized states, namely, $`\left|\psi |\psi ^{}\right|^2`$. It is directly related with the natural distance in the space of rays, that is, the angle between rays. The geometrical meaning of this distance is best perceived by considering the complex projective space of rays, where it is called the Fubini-Study metric .111It is very interesting to realize that the distance between the distribution probabilities defined by the outcomes of possible measurements of observables that distinguish the states $`\psi `$ and $`\psi ^{}`$ coincides with the statistical distance in the classical sense of Fisher . Schumacher maximization of the fidelity uses the Schmidt decomposition of the entangled state .

### 2.3 Renormalization group and information theory

The problem of transforming one quantum state into another while preserving its fundamental features is reminiscent of an operation performed with the quantum renormalization group (in this connection, see Ref. ). This operation does not intend to reach near-perfect fidelity but it is desirable that it reach as much fidelity as possible. We shall see that a particular formulation of the renormalization group, namely, the density matrix renormalization group, comes close to applying the concepts of quantum information theory, and in a similar way to Rindler space quantization . In fact, the construction of the density matrix renormalization group algorithm is based on the Schmidt decomposition, in parallel to Schumacher fidelity maximization.

Let us consider an entangled quantum state $`|\psi `$ for which we seek an optimal reduced representation $`|\stackrel{~}{\psi }`$ (a sort of quantum coding). That is, we must find a projection

$$|\psi =\underset{aA}{}\psi _{aA}|a|A|\stackrel{~}{\psi }=\underset{iA}{}\psi _{iA}|i|A$$
(9)
to a subspace spanned by reduced basis $`\{i\}`$ such that the distance $`S=||\stackrel{~}{\psi }|\psi |^2`$ is minimized. This is the problem that S. White met in looking for an improved “real space” quantum renormalization group algorithm: In a renormalization group step, one must reduce the number of block states, but in a way that accounts for the influence of the rest of the system (that is, with no introduction of arbitrary boundary conditions) . This amounts to consider the entanglement of the block and to formulate the above described problem. White found the solution in terms of the singular value decomposition, a well-known numerical algorithm, equivalent to the Schmidt decomposition of the entangled state and, hence, leading to discarding the smallest eigenvalues of the block density matrix; in other words, $`|i`$ are the eigenstates of $`\rho `$ with largest eigenvalues. The proper formulation of the density matrix renormalization group in terms of quantum information concepts has been provided recently .

#### 2.3.1 Density matrix renormalization group algorithm

Method for a $`1D`$ quantum system on a chain (e.g., a chain of oscillators):

1. Select a sufficiently small “soluble” block $`[0,L]`$:
2. Reflect the block on the origin:
3. Compute the ground state.
4. Compute the density matrix of the block $`[0,L]`$.
5. Discard eigenstates with smallest eigenvalues.
6. Add one site next to the origin.
7. Go to 2.

One has to adjust this procedure in such a way that the iteration keeps the Hilbert space size approximately constant. The procedure can be performed algebraically for a chain of coupled harmonic oscillators . Otherwise, it has to be performed numerically. In the continuum limit, the connection with the scalar field theory half-space density matrix previously calculated is clear. In particular, the property of “angular waves” of having vanishing wave-length at $`x=0`$ pointed out above provides an explanation of the form in which the density matrix renormalization group solves the boundary condition problem, that is, because of the concentration of quantum states near the boundary. To be precise, using the eigenfunctions of $``$ instead of free waves, we have a basis in which the region close to $`x=0`$ (the boundary point) is more accurately represented than the region far from it when we cut off the higher $`\mathrm{}`$ eigenfunctions.

## 3 Generalizations and applications

### 3.1 Geometric entropy

We have seen that the half line density matrix of a field theory has a geometric interpretation in Rindler space. Furthermore, the entropy of black holes can be understood as a generalization to a more complicated geometry. We may wonder if further generalizations are possible.

In this connection, we recall the notion of “geometric entropy”, introduced by C. Callan and F. Wilczek , as the entropy “associated with a pure \[global\] state and a geometrical region by forming the pure state density matrix, tracing over the field variables inside the region to create an ‘impure’ density matrix”. Of course, their motivation was the earlier suggestion that back-hole entropy is of quantum-mechanical entanglement origin . They proposed a generalization to different topologies but, actually, they only computed the Rindler space case, discussing the divergence of the entropy at the horizon  (the UV divergence of this type of entropy had been discussed in general in Ref. ).

A different notion of geometric entropy can be deduced by purely geometrical means from the presence of horizons, namely, as associated with a spacetime topology that does not admit a trivial Hamiltonian foliation . This type of topology prevents unitary evolution and leads to mixed states.

In fact, it is only the second type of entropy that embodies the famous “one-quarter area law” for black holes, due to its origin in purely relativistic concepts (this was demonstrated for an earlier relativistic notion of entropy in Ref. ). On the contrary, the quantum notion of geometric entropy for a field theory involves UV divergences and needs renormalization before a comparison with the relativistic notion can be made (see the discussion in Ref. ).

### 3.2 Application to cosmology

The generalization of the concept of black-hole entropy to the de Sitter space and, hence, to cosmology is relatively old . Of course, the concept of entropy in its traditional thermodynamical sense has been crucial in explaining the dynamics of inflation (now a standard paradigm), namely, in accounting for the reheating process (entropy generation). However, the relations between the traditional view and the one associated to quantum entanglement may lead to further insight, when they are properly formulated in the context of quantum information theory.

In particular, we can regard the generation of entropy and fluctuations in de Sitter spacetime as a fundamentally quantum process leading to the celebrated Harrison-Zeldovich scale invariant spectrum of Gaussian fluctuations. If we consider the initial state for inflation as a pure quantum state (that some theory of quantum gravity will hopefully characterize some day), the de Sitter space horizon induces decoherence of the modes which cross it, irrespective of the actual inflaton dynamics and, therefore, of the details of the reheating process. This decoherence consists of a randomization of the phases of the present quantum fields and naturally produces thermal Gaussian fluctuations. Moreover, the symmetry of de Sitter space implies that each mode has the same physical size as it crosses the horizon, leading to the Harrison-Zeldovich power spectrum:

$$P(k)\left|\left(\frac{\delta \rho }{\rho _0}\right)_k\right|^2k.$$

Actually, this spectrum is given by the equipartition theorem at the corresponding Hawking temperature $`T=H/(2\pi )`$ ($`H=\mathrm{\Lambda }/3`$) .

In this reasoning, one could also consider the initial state to be mixed (e.g., a thermal state) and apply the argument for entropy growth exposed in subsection 1.1.

One may wonder if a density-matrix type renormalization group could help with the dynamics. As long as the fluctuations are Gaussian, we are in a similar situation to the case of harmonic oscillators commented above, which makes any renormalization group superfluous. However, as is well known, the gravitational instability produces non-linear evolution and leads to phase correlations. Therefore, the ideas presented here may be useful in setting up a renormalization group for the study of non-Gaussian fluctuations and its non-linear evolution. In particular, the Wilson or exact renormalization group irreversibility properties and can be connected with other methods of analysis of the non-linear evolution (for preliminary steps in that direction, see Ref. ).

Acknowlegments

I am grateful to Juan Poyatos for bringing some relevant references to my attention.


text:
# Spin-dependent tunnelling through a symmetric barrier

## Abstract

The problem of electron tunnelling through a symmetric semiconductor barrier based on zinc-blende-structure material is studied. The $`k^3`$ Dresselhaus terms in the effective Hamiltonian of bulk semiconductor of the barrier are shown to result in a dependence of the tunnelling transmission on the spin orientation. The difference of the transmission probabilities for opposite spin orientations can achieve several percents for the reasonable width of the barriers.

Lately spin polarized electron transport in semiconductors attracts a great attention. spintronic  One of the major problems of general interest is a possibility and methods of spin injection into semiconductors. A natural way to achieve spin orientation in experiment is the injection of spin polarized carriers from magnetic materials. Although significant progress has been made recently, Zhu ; Hammar ; Hanbicki ; Safarov  reliable spin-injection into low-dimensional electrons systems is still a challenge. Schmidt et al. pointed out that a fundamental obstacle for electrical injection from ferromagnetic into semiconductor was the conductivity mismatch of the metal and the semiconductor structure. Schmidt  However, Rashba showed that this problem could be resolved by using tunnelling contact at the metal-semiconductor interface. Rashba2  On the other hand Voskoboynikov et al. Vosk1  proposed that asymmetric non-magnetic semiconductor barrier itself could serve as a spin filter. It was demonstrated that spin-dependent electron reflection by inequivalent interfaces resulted in the dependence of the tunnelling transmission probability on the orientation of electron spin. This effect is caused by interface-induced Rashba spin-orbit coupling Rashba  and can be substantial for resonant tunnelling through asymmetric double-barrier Vosk2 ; Silva  or triple-barrier Koga  heterostructures. However, in the case of symmetric potential barriers, the interface spin-orbit coupling does not lead to a dependence of tunnelling on the spin orientation.

In this communication we will show that the process of tunnelling is spin dependent itself. We demonstrate that a considerable spin polarization can be expected at tunnelling of electrons even through a single symmetric barrier if only the barrier material lacks a center of inversion like zinc-blende structure semiconductors. The microscopic origin of the effect is the Dresselhaus $`k^3`$ terms Dresselhaus  in the effective Hamiltonian of the bulk semiconductor of the barrier.

We consider the transmission of electrons with the initial wave vector $`𝒌=(𝒌_{},k_z)`$ through a flat potential barrier of height $`V`$ grown along $`z[001]`$ direction (see Fig. 1); $`𝒌_{}`$ is the wave vector in the plane of the barrier, and $`k_z`$ is the wave vector normal to the barrier pointing in the direction of tunnelling. The electron Hamiltonian of the barrier in effective mass approximation contains the spin-dependent $`k^3`$ term (Dresselhaus term) Dresselhaus

$$\widehat{H}_D=\gamma \left[\widehat{\sigma }_xk_x(k_y^2k_z^2)+\widehat{\sigma }_yk_y(k_z^2k_x^2)+\widehat{\sigma }_zk_z(k_x^2k_y^2)\right],$$
(1)
where $`\widehat{\sigma }_\alpha `$ are the Pauli matrices, $`\gamma `$ is a material constant (see Table 1), and the coordinate axis $`x,y,z`$ are assumed to be parallel to the cubic crystallographic axis $`[100]`$, $`[010]`$, $`[001]`$, respectively. In the case of tunnelling along $`z`$ one should consider $`k_z`$ in the Hamiltonian as an operator $`i/z`$. We assume the kinetic energy of electrons to be substantial smaller than the barrier high $`V`$, then the Hamiltonian (1) is simplified to

$$\widehat{H}_D=\gamma (\widehat{\sigma }_xk_x\widehat{\sigma }_yk_y)\frac{^2}{z^2}.$$
(2)
One can note that essentially $`\widehat{H}_D`$ induces a spin-dependent correction to the effective electron mass along $`z`$ in the barrier.

|  | GaSb | InAs | GaAs | InP | InSb |
| --- | --- | --- | --- | --- | --- |
| $`\gamma `$, eV$``$Å3 | 187 | 130 | 24 | 8 | 220 |
| $`m^{}/m_0`$ | 0.041 | 0.023 | 0.067 | 0.081 | 0.013 |

Table 1: Parameters of band structure of various A3B5 semiconductors. PT ; IP

The Hamiltonian (2) is diagonalized by spinors

$$\chi _\pm =\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ \text{e}^{i\phi }\end{array}\right),$$
(3)
which correspond to the electron states ”$`+`$” and ”$``$” of the opposite spin directions. Here $`\phi `$ is the polar angle of the wave vector $`𝒌`$ in the plane $`xy`$, being

$$𝒌=(k_{}\mathrm{cos}\phi ,k_{}\mathrm{sin}\phi ,k_z).$$
(4)
Transmission probabilities for the electrons of eigen spin states ”$`+`$” and ”$``$” (3) are different due to spin-orbit term (2). The orientations of spins $`𝒔_\pm `$ in the states ”$`+`$” and ”$``$” depend on the in-plane wave vector of electrons and are given by

$$𝒔_\pm =(\mathrm{cos}\phi ,\pm \mathrm{sin}\phi ,\mathrm{\hspace{0.25em}0}).$$
(5)
Fig. 2 demonstrates the orientations of spins $`𝒔_+`$ and $`𝒔_{}`$ for various directions of the in-plane electron wave vector $`𝒌_{}`$, i.e. as a function of polar angle $`\phi `$. If $`𝒌_{}`$ is directed along a cubic crystal axis ($`[100]`$ or $`[010]`$) then the spins are parallel (or antiparallel) to $`𝒌_{}`$, while $`𝒔_\pm `$ are perpendicular to $`𝒌_{}`$ if the in-plane wave vector is directed along the axis $`[1\overline{1}0]`$ or $`[110]`$.

Electrons of the eigen spin states ”$`+`$” and ”$``$” propagate through the barrier with conserving of the spin orientation. Since the wave vector in the plane of the barrier $`𝒌_{}`$ is fixed, wave functions of the electrons can be written in the form

$$\mathrm{\Psi }_\pm (𝒓)=\chi _\pm u_\pm (z)\mathrm{exp}(i𝒌_{}𝝆),$$
(6)
where $`𝝆=(x,y)`$ is a coordinate in the plane of the barrier. The function $`u(z)`$ in the regions I (incoming and reflected waves, see Fig. 1), II and III (transmitted wave) has the form

$$u_\pm ^{(I)}(𝒓)=\left[\mathrm{exp}(ik_zz)+r_\pm \mathrm{exp}(ik_zz)\right],$$
(7)
$$u_\pm ^{(II)}(𝒓)=\left[A_\pm \mathrm{exp}(q_\pm z)+B_\pm \mathrm{exp}(q_\pm z)\right],$$

$$u_\pm ^{(III)}(𝒓)=t_\pm \mathrm{exp}(ik_zz),$$

respectively. Here $`t_\pm `$ and $`r_\pm `$ are the transmission and reflection coefficients for spin states $`\chi _\pm `$, respectively, and the wave vectors under the barrier $`q_\pm `$ are given by

$$q_\pm =q_0\left(1\pm \gamma \frac{2m_2k_{}}{\mathrm{}^2}\right)^{1/2},$$
(8)
where $`q_0`$ is the reciprocal length of decay of the wave function in the barrier for the case when the spin-orbit interaction (1) is neglected

$$q_0=\sqrt{\frac{2m_2V}{\mathrm{}^2}k_z^2\frac{m_2}{m_1}k_{}^2\left(\frac{m_2}{m_1}1\right)},$$
(9)
and $`m_i`$ ($`i=1,2`$) are the effective masses outside and inside the tunnelling barrier, respectively. Taking into account the boundary conditions, which require that

$$u_\pm \text{and}\frac{1}{m}\frac{u_\pm }{z}$$
(10)
are continuous at the interfaces, a system of linear equations for $`t_\pm `$, $`r_\pm `$, $`A_\pm `$, and $`B_\pm `$ can be derived. Solution of the system allows one to calculate the coefficients of the transmission $`t_\pm `$. For the real case $`m_2k_{}\gamma /\mathrm{}^21`$ they are derived to be

$$t_\pm =t_0\mathrm{exp}\left(\pm \gamma \frac{m_2k_{}}{\mathrm{}^2}aq_0\right),$$
(11)
where $`t_0`$ is the transmission coefficient when the spin-orbit interaction (1) is neglected,

$$t_0=4i\frac{m_2}{m_1}\frac{k_zq_0}{\left(q_0ik_zm_2/m_1\right)^2}\mathrm{exp}(aq_0),$$
(12)
$`a`$ is the width of the barrier. The general problem of tunnelling of an electron with arbitrary initial spinor $`\chi `$ can be solved by expanding $`\chi `$ to the eigen spinors $`\chi _\pm `$.

It is convenient to introduce a polarization efficiency $`𝒫`$ that determines the difference of tunnelling transmission probabilities for the spin states ”$`+`$” and ”$``$” through the barrier

$$𝒫=\frac{|t_+|^2|t_{}|^2}{|t_+|^2+|t_{}|^2}.$$
(13)
In our case it has the form

$$𝒫=\mathrm{tanh}\left(2\gamma \frac{m_2k_{}}{\mathrm{}^2}aq_0\right).$$
(14)
At a given initial wave vector of electrons, $`𝒌`$, the polarization efficiency drastically increases with the strength of the Dresselhaus spin-orbit coupling $`\gamma `$ (see Eq. 14), and the barrier width $`a`$. However while increase the barrier width $`a`$ increases the tunnelling efficiency one should keep in mind that the barrier transparency decreases simultaneously (see Eq. (12)). In Fig. 3 the efficiency $`𝒫`$ and the barrier transparency $`|t_0|^2`$ are plotted as a function of the barrier width, $`q_0a`$, for various barrier materials.

The material parameters $`\gamma `$ and effective mass $`m^{}=m_2`$ used in the calculations are given in the Table 1. One can see that for a reasonable set of parameters given in the figure caption it is possible to achieve spin polarizations of several percents. Tunnelling barriers prepared on the basis of GaSb or its solutions seem to be the most efficient barrier materials for spin selective tunnelling because of the large value of the product $`\gamma m^{}`$.

The polarization strongly depends on the electron wave vector $`k_{}`$ parallel to the barrier (see Eq. (14)). This result suggests a device for spin injection into quantum wells. Let’s assume two quantum wells separated by a tunneling barrier, and a current flowing along one of the quantum wells. The in-plane current results in non-zero average electron wave vector $`k_{}`$ and, due to the considered effect, in a spin polarization of carriers.

In conclusion, we have demonstrated that the $`k^3`$ Dresselhaus terms in the effective Hamiltonian of semiconductors lacking a center of inversion yield a considerable spin polarization of electrons tunnelling through barriers. The effect could be employed for creating spin filters, eg. on the base of type-II strained heterostructures like InAs/GaSb, InSb/GaSb and GaSb/GaAs.

The work was supported by the RFBR, the Presidium of the RAS, and grants of the DFG and INTAS.