format: text

text:
# Solution stability and variability in a simple model of globular proteins

## I Introduction

Proteins are linear heteropolymers: they are linear sequences of monomers, each of which is one of twenty different types. Different proteins have different sequences of amino acids. These differences allow proteins to perform the huge range of tasks they do in living cells. But this does not mean that 2 proteins that do the same job necessarily have the same sequence. For example, many organisms have enzymes called adenylate kinases which perform essentially the same job in the cytoplasm of each organism. But the amino acid sequences of adenylate kinases vary very widely, even though they are all doing the same job in more-or-less the same milieu. Below are the amino-acid sequences of the adenylate kinases of two prokaryotes proteins . First that of Escherichia coli

$$\begin{array}{c}\text{MRIILLGAPGAGKGTQAQFIMEKYGIPQISTGDMLRAAVKSGSELGKQAK}\hfill \\ \text{DIMDAGKLVTDELVIALVKERIAQEDCRNGFLLDGFPRTIPQADAMKEAG}\hfill \\ \text{INVDYVLEFDVPDELIVDRIVGRRVHAPSGRVYHVKFNPPKVEGKDDVTG}\hfill \\ \text{EELTTRKDDQEETVRKRLVEYHQMTAPLIGYYSKEAEAGNTKYAKVDGTK}\hfill \\ \text{PVAEVRADLEKILG}\hfill \end{array}$$

and secondly that of Vibrio cholerae

$$\begin{array}{c}\text{MRIILLGAPGAGKGTQAQFIMEKFGIPQISTGDMLRAAIKAGTELGKQAK}\hfill \\ \text{AVIDAGQLVSDDIILGLIKERIAQADCEKGFLLDGFPRTIPQADGLKEMG}\hfill \\ \text{INVDYVIEFDVADDVIVERMAGRRAHLPSGRTYHVVYNPPKVEGKDDVTG}\hfill \\ \text{EDLVIREDDKEETVRARLNVYHTQTAPLIEYYGKEAAAGKTQYLKFDGTK}\hfill \\ \text{QVSEVSADIAKALA}\hfill \end{array}$$

where the sequences are given as a sequence of the 1-letter codes for the amino acids of which they are made. The first amino acid is an M (Methionine), the second is an R (Arginine) and so on. The sequence is read as English text, from top left to bottom right. See any molecular biology or biochemistry textbook thecell ; stryer ; voet  for an introduction to amino acids and proteins. Note that there are many differences between the sequences! The amino-acid sequences of proteins are very different while keeping the function. Also, we picked adenylate kinases only in order to have a concrete example, it is a general property of proteins. The function of adenylate kinases is irrelevant to our discussion of stability, beyond the fact that they function as enzymes as monomers in solution inside cells. Here we will concentrate entirely on globular proteins, the proteins that exist in solution not embedded in membranes.

Now, the simplest thing to do when faced with this radical difference in sequence without a corresponding difference in function is to ignore it. To assume that the 2 proteins interact and behave in a very similar manner. But do they? As they both function as proteins inside the cytoplasm of bacteria they are both clearly soluble and do not stick to things they should not stick to in vivo. However, this does not mean that their solubilities, for example, are necessarily equal. Both their solubilities are sufficient to allow them to function but one may exceed the minimum solubility by a large margin and one by a small margin. It would be of interest to know what these margins are and how they vary from protein to protein, not only because we wish to understand how proteins function and have evolved in vivo, but to help us process, purify, and crystallise proteins. If a protein is only marginally soluble in the conditions in vivo then it may aggregate when its environment (salt concentrations, temperature etc.) are altered. We would like to understand and to be able to predict, the variability of properties, such as solubility, of proteins.

We will focus on the stability of solutions of proteins in their native state, i.e., we assume that the protein has folded into its native state and remains there. Thus we consider only folded proteins sticking together due to their surfaces attracting each other, not proteins partially unfolding and then aggregating due to the hydrophobic regions of the protein exposed by unfolding, attracting each other. So, our proteins will always be compact objects, more like colloidal particles than conventional polymers. This allows us to avoid the complex problem of protein (un)folding. Effectively, we assume that proteins such as the adenylate kinases of E. coli and V. cholerae differ only in their surfaces. Replacing one surface amino acid in the chain by another then changes only the surface and through that the protein-protein interaction. If a hydrophobic amino acid replaces a hydrophilic amino acid in a position on the chain where the chain is at the protein’s surface, then we expect the surface to become more sticky, which would tend to decrease the second virial coefficient, whereas replacing a hydrophobic amino acid by a hydrophilic one should have the opposite effect. For simplicity, instead of having 20 different types of amino acids at the surface, we use a model whose surface is described by bits which have only 2 values: hydrophobic and hydrophilic. This is a rather gross approximation, the amino acids vary widely in size, some are charged, but we want the simplest possible model. The model is an extension of that considered in Ref. sear02, . A protein molecule is modeled by a cube, whose 6 faces interact with a short-ranged attraction, which is here determined by a sequence of $`n_B`$ bits. In Ref. sear02,  the interaction between faces was taken to be a random variable; we will discuss the differences between that model and the more complex one considered here in the conclusion. A schematic of the model is shown in Fig. 1.

We have talked of our model proteins being soluble in vivo. Real proteins have evolved to be so. The cytoplasm of bacteria such as E. coli and V. cholerae is very complex: bacteria typically have a few thousand different proteins vibnote , and any one of these proteins is then surrounded by thousands of different proteins, as well as RNA and DNA, small molecules such as nucleotides etc.. An individual enzyme must be soluble in the sense that it does not stick too strongly to not only other proteins of the same type but those of all the other types, as well as not binding to the RNA, DNA, etc.. In future work, we will address this problem, but here we will keep things simple and consider only interactions between model proteins of one type. We will calculate the second virial coefficient only for the interaction of 2 model protein molecules of the same type. This is not realistic for a enzyme in a bacterial cell as an individual enzyme will be present at rather low concentrations, even though the total protein concentration in bacterial cells is around 20% by volume neidhardt ; thecell . It is however, a good place to start, and is realistic for a few exceptional cells, such as our red blood cells which contain very high concentrations of a single protein: hemoglobin. Future work will address this issue and will also look at proteins which bind to other proteins, as many proteins do.

The proteins whose sequences we gave in the first paragraph are presumably orthologs: they are both descended from a common ancestral protein but have evolved independently, keeping their function the same, since the E. coli and V. cholerae lineages separated. The fact that proteins with the same function, but that have evolved independently in different species, can have very different amino-acid sequences, is well known. Paralogs, proteins created by duplication of a gene, also start with identical sequences but have sequences that diverge with time. The differences are believed to have arisen via random mutations which are not rejected by natural selection because they are not actively deleterious (to the survival of the organism) but also do not have any selective advantage. This theory of mutations changing the amino-acid sequences of proteins without improving or reducing its ability to function is called the theory of neutral evolution kimura68 ; king69 ; ayala97 . The constraints placed on this neutral evolution by the requirement on the protein to fold have been considered bastolla00 ; bastolla02 ; deeds03 ; tang00 ; bornberg99 ; taverna02 , but not those due to the requirement of the protein to be soluble. The constraints placed on the sequences of RNA by the requirement to be functional and the evolution of these sequences, are analogous to the constraints on the sequences of, and evolution of proteins. They have been extensively studied and in many respects are rather better understood, essentially because RNA is simpler than protein. See the review of Higgshiggs . However, there has been some work which has considered protein-to-protein variability sear02 ; sear03 , see also Ref. rosenbaum99,

We will generate our model proteins at random (subject to the solubility constraint) and assume that neutral evolution of proteins is close to a random walk from one sequence to another. This random walk occurs in what is often called ‘protein space’ smith70 , with each sequence a unique point in this space and 2 sequences neighbours if 1 of them can be transformed into the other by a single mutation. This protein space is vast. The set of soluble proteins exists in this protein space as a set of points, 1 for each soluble protein. A schematic of the protein spaces of proteins and model proteins is shown in Fig. 2. It is only very schematic, the space is huge and many dimensional. In each case the arrows represent a single mutation changing a protein into a neighbouring protein. Below, we will generate random walks for our model proteins, and these will sample all soluble states with equal probability. When we come to applying our results to real, not model, proteins, we will have to assume that neutral evolution also samples proteins which are soluble with reasonably uniform probabilities.

In the following section, we will perform a simple analysis of sequence data, to look at variations in the number of hydrophobic amino acids. The model is defined in section III, and the stability of its solutions estimated and discussed in section IV. The last section is a conclusion.

## II Analysis of sequence data

The sequences of the adenylate kinases of E. coli and V. cholerae are both of viable enzymes, they are soluble in vivo and catalyse a reaction. Looking at them, an obvious question to ask is: How many sequences of amino acids are there, that fold up to form viable adenylate kinases? Both adenylate kinases have 214 amino acids. As there are 20 types of amino acids there are $`20^{214}10^{278}`$ different amino acid sequences of 214 amino acids. An enormous number, of which presumably the vast majority do not fold into a unique native state, let alone are soluble and act as a catalyst. But it seems likely that the number of possible amino acid sequences that correspond to viable adenylate kinases is huge.

A database at SWISSPROT swissprot ; spref , called PROSITE prosite ; prositeref , identified 152 amino-acid sequences as belonging to the adenylate-kinase family of proteins (PROSITE accession number PS00113). It did so by locating the amino acids of the active site of an adenylate kinase prosite ; prositeref . 104 of these sequences are from prokaryotes, of which we eliminate 4 sequences as they contain less than 100 amino acids and are presumably not complete proteins. This leaves 100 adenylate kinases; 2 of these kinases are the ones whose sequences are in the first paragraph. We can calculate the fraction of the amino acids of these adenylate kinases that are hydrophobic, $`h`$, and plot this against the length of the sequence, $`M`$: the total number of amino acids in the sequence. The results are shown as a scatter plot, Fig. 3. The 9 amino acids G, A, V, L, I, M, P, F and W, are taken to be hydrophobic, and the remaining 11 to be hydrophilic. Here each amino acid is represented by its 1-letter code: G for glycine, A for alanine, etc.. The 9 hydrophobic amino acids are those whose side chains are classified as nonpolar in Ref. voet,  (Table 4-1, p58). There is some arbitrariness in where the dividing line is drawn between hydrophobic and hydrophilic amino acids, but different dividing lines give rather similar spreads in $`h`$.

For the present work, the key observation is that the fraction of an adenylate kinase’s amino acids which are hydrophobic varies from protein to protein, as do other properties such as their net charge sear03 ; unpub . In section IV we will find that for our model, with a constraint imposed that model proteins are soluble, there is scatter in the fraction of its bits that are hydrophobic.

## III Model

The model is chosen to be as simple and as generic as possible, while having interactions which are mediated by surface patches whose interactions are a function of sequences or string of bits. The protein-protein interactions then depend on the values of these bits, some sets of values give proteins which strongly attract each other while other sets give proteins which largely repel each other. This is perhaps the simplest model of a globular protein which allows for mutations. Within the model these mutations flip one of the bits, a model of a mutation which converts a surface residue from a hydrophobic amino acid to a hydrophilic amino acid, or vice versa. A schematic of the model is shown in Fig. 1. An amino acid of a protein is called a residue.

The model protein is a cube, with each of its 6 faces having a single patch sear02 . The lattice is cubic and each protein occupies 8 lattice sites arranged 2 by 2 by 2, see Fig. 1. We make the model 2 sites across to reduce the range of the attraction, which is 1 site, to half the diameter of the hard core. The model proteins can rotate, and so have 24 distinct orientations. Each of the 6 faces of the cube has a patch, labeled $`i=1`$ to 6, with patches 1 to 4 clockwise around a loop of 4 of the faces, and patches 5 and 6 on the remaining 2 faces. The interactions between model proteins are pairwise additive and consist of 2 parts. The first is simply an excluded-volume interaction: 2 proteins cannot overlap. The second is that if the faces of 2 proteins are in contact there is an energy of interaction between the 2 touching patches of the 2 proteins. By in contact we mean that the faces must overlap completely otherwise the energy of interaction is taken to be zero. Also, the model is such that the energy of interaction between two touching patches is a constant which does not change when the two proteins are rotated about the axis joining their centres. The touching patches are those on the faces of the 2 proteins that face each other. This is all as in Ref. sear02, , the difference is in how the interaction energy of a pair of patches $`i`$ and $`j`$, $`u_{ij}`$, is specified.

How a patch interacts is specified by a sequence or string of $`n_B`$ bits. If a bit has a value of 1 then the bit is said to be hydrophilic or polar, whereas if it has a value of 0 then it is hydrophobic. The interaction energy of a pair of touching patches, $`i`$ and $`j`$, is then given by

$$u_{ij}=ϵ\underset{\alpha =1}{\overset{n_B}{}}\left(b_\alpha ^{(i)}1\right)\left(b_{1+n_B\alpha }^{(j)}1\right),$$
(1)
where $`b_\alpha ^{(i)}`$ is bit number $`\alpha `$ of patch $`i`$. $`ϵ`$ is the interaction energy of 2 hydrophobic bits. We use energy units such that the thermal energy $`k_BT=1`$. Thus to calculate the interaction the string of bits of 1 of the patches is reversed and then the energy is just the sum of the number of pairs of corresponding bits where both bits are 0, are hydrophobic. The only interaction is between 2 hydrophobic bits; there is no hydrophobic-hydrophilic or hydrophilic-hydrophilic interaction. The reason one of the strings is reversed is that if this is not done then the interaction between like patches, $`j=i`$, is just $`ϵ`$ times the number of 0s in $`i`$’s string. Reversing the strings removes this problem in a simple way. Of course, the interactions form a symmetric square matrix, $`u_{ij}=u_{ji}`$. Each of the 6 patches is taken to be labeled and so distinguishable, i.e., we take a pair of proteins where one protein can be obtained from the other by swapping a pair of the strings of bits, as 2 different proteins.

Thus, a protein is specified by giving values to the $`6`$ strings of $`n_B`$ bits, and so there are $`2^{6n_B}`$ possible different proteins. For all but rather small values of $`n_B`$, this is a very large number of possible proteins, e.g., for $`n_B=18`$, we have $`3\times 10^{32}`$ different model proteins. This is however, much smaller than the number of possible real proteins. Most of the calculations have been done for $`n_B=18`$, with a few for $`n_B=12`$, for comparison. We choose $`n_B=18`$ as being a sensible number as then the total number of bits which describe the surface is 108. Adenylate kinases, for example, have around 200 amino acids, of which about half are on the surface. Thus, we have about 1 bit per surface amino acid. Our model proteins can be thought of as existing in ‘protein space’ with each possible protein represented by a point in this space, and each protein has $`6n_B`$ neighbours, each of which is obtained by flipping 1 of the bits of the protein, see Fig. 2.

The second virial coefficient $`B_2`$ of our lattice model is given by sear02

$$B_2=\frac{1}{2}\left[27\frac{1}{6}\underset{i=1}{\overset{6}{}}\underset{j=1}{\overset{6}{}}\left(\mathrm{exp}\left(u_{ij}\right)1\right)\right],$$
(2)
where the first term inside the brackets comes from excluded-volume interactions and the second from the interactions between touching patches. The number 27 comes from the fact that each model protein excludes other proteins from a cube of 3 by 3 by 3 lattice sites. Thus, in the high temperature limit $`B_2=B_{2hc}=27/2`$. The sums over 24 orientations reduce to sums over 6 orientations as rotating either of the 2 molecules around the axis joining their centres does not change the energy. The factor in front of the double sum is a normalisation factor of $`1/36`$ times the 6 possible lattice sites that one molecule can occupy and be adjacent to the other molecule.

## IV Stability of solutions

Unless $`ϵ`$ is small, many of the $`2^{6n_B}`$ model proteins strongly attract each other leading to condensation, gelation, and possibly crystallisation. By condensation we mean the formation of coexisting dilute and concentrated protein solutions, as have been studied extensively for the protein lysozyme muschol97 ; vliegenthart00 ; rosenbaum99 . Only a fraction of the model proteins are viable in the sense that they are stable as single phase solutions. Clearly proteins cannot condense in vivo without severely impairing the organism’s function.

The attractions affect the phase behaviour through and can be measured by, the second-virial coefficient. In the absence of attractions the second-virial coefficient is approximately 4 times the volume of a particle (assuming the particle is not too anisotropic). Attractions decrease its value until eventually the pressure does not increase monotonically but decreases over a range of densities due to the negative virial coefficient; a van der Waals loop forms. If we impose the constraint that the second virial coefficient be above a certain value, where we believe the pressure will be a monotonic function of density, we can quantify what fraction of our model proteins satisfy this constraint and so have solutions which are stable.

We insist that the reduced second virial coefficient satisfy $`B_2/B_{2hc}1`$, in order for the protein to be viable. The fraction of proteins which are viable, according to this criterion, is denoted by $`f_v`$. It is determined by generating proteins at random, setting each bit to be 0 or 1 with equal probabilities, and finding the fraction with $`B_2/B_{2hc}1`$. See the appendix for further details of the computations. The value of $`B_2/B_{2hc}`$ at the critical point, the highest point on the curve separating the 1 and 2-phase regions of a phase transition into coexisting solutions, is typically a little less than $`1`$, unless the attraction is very anisotropic. For the canonical model, hard spheres plus a long-range attraction, the critical point occurs when $`B_2/B_{2hc}=1.65`$, and provided the attraction remains isotropic this value changes little even if the attraction is made quite short ranged vliegenthart00 . If the attraction is very anisotropic then $`B_2/B_{2hc}`$ can (depending a little on the precise nature of the anisotropy) be much more negative at the critical point sear99 ; curtis01 ; kern03 , but for simplicity we insist on $`B_2/B_{2hc}`$ being above a fixed value for all our proteins, regardless of how anisotropic their attractions are. Crystallisation out of not-too-concentrated solutions also requires as a minimum, attractions of about the strength required to make $`B_2/B_{2hc}`$ around $`1`$. The propensity to crystallise depends on the details of the attraction, for work on the earlier version of this model with random values of the patch-patch attractions, see Ref. sear02, .

Results are shown, as a function of $`ϵ`$, for $`n_B=12`$ and 18, in Fig. 4. As might have been expected, as $`ϵ`$ increases, the fraction of viable proteins decreases exponentially, but note that even for $`n_B=12`$ and $`ϵ=2`$, there are still $`7.1\times 10^9`$ viable proteins, a very large number. Partly, what is happening is that as $`ϵ`$ increases then fewer and fewer hydrophobic bits are allowed, and as the fraction of bits that are hydrophobic decreases, then the number of possible proteins decreases: there are many possible proteins with close to half their bits 0s and half 1s, but only one with all its bits equal to 1. Partly, what happens is that correlations are introduced between the hydrophobic bits in the strings. The hydrophobic bits tend to avoid each other, e.g., if all 6 strings have all their bits from 1 to $`n_B/2`$ (assuming $`n_B`$ is even) hydrophilic, then any or all of their bits from $`n_B/2+1`$ to $`n_B`$ may be hydrophobic without there being any attractions. Thus here the hydrophobic bits avoid each other, in order to avoid the attractive interactions which make the second virial coefficient negative and thus violate our solubility condition.

We can measure both these effects by defining 2 quantities. The first is the mean fraction of bits which are 0, are hydrophobic. Denoting this by $`h`$, it is defined by

$$h=\frac{1}{6n_B}\underset{i=1}{\overset{6}{}}\underset{\alpha =1}{\overset{n_B}{}}\left(1b_\alpha ^{(i)}\right).$$
(3)
The average denoted by $``$ is over proteins which satisfy our criterion for the stability of the solution. We use $`h`$ to denote both the fraction of bits in our model proteins that are hydrophobic, and the fraction of residues in real proteins that are hydrophobic. A measure of the correlation between the probability that a bit $`\alpha `$ is hydrophobic, and that the bit $`1+n_B\alpha `$ with which it interacts is also hydrophobic is denoted by $`hh_p`$, and is defined by

$$hh_p=\frac{1}{36n_Bh^2}\underset{i=1}{\overset{6}{}}\underset{j=1}{\overset{6}{}}\underset{\alpha =1}{\overset{n_B}{}}\left(1b_\alpha ^{(i)}\right)\left(1b_{1+n_B\alpha }^{(j)}\right).$$
(4)
We have plotted both quantities in Fig. 5. The model has $`n_B=18`$ bits and the quantities are plotted as a function of $`ϵ`$. As $`ϵ`$ increases, the fraction of bits which can be hydrophobic without the second virial coefficient becoming too negative decreases. Also, the anticorrelations between a bit being hydrophobic and the bit with which it interacts being also hydrophobic increases. If there were no correlation between the states of the 2 bits then $`hh_p=1`$, which is true for $`ϵ=0`$, but this function decreases as $`ϵ`$ increases. If a bit is hydrophobic the bit with which it interacts is less likely to be hydrophobic. We have shown results just for $`n_B=18`$ but results for other numbers of bits are similar.

We only constrain the second virial coefficient to be above a certain value, we do not constrain its precise value. As the second virial coefficient is a function of the number of hydrophobic bits on its 6 faces and as this number is an integer between 0 and $`n_B`$, the second virial coefficient can only take one of a set of values, and so the probability density function for $`B_2/B_{2hc}`$ is a set of delta functions. We have plotted these as spikes, with the height of each spike set to the probability that $`B_2/B_{2hc}`$ has this value. We can see that the most likely values of the reduced second virial coefficient are near the minimum allowed value of $`1`$. This is simply because there are many more sets of strings with close to half the bits hydrophobic than there are with most of the bits hydrophilic, and the proteins with close to half the bits hydrophobic have very large and negative second virial coefficients. There is only one protein with all 108 bits hydrophilic but the number of proteins which have 9 hydrophobic and 9 hydrophilic bits on each face is $`(18!/9!^2)^610^{28}`$. The probability distribution function for all possible proteins (including those with $`B_2/B_{2hc}<1`$) is sharply peaked at a value much less than 1, for $`n_B=18`$ and $`ϵ=2`$, and Fig. 6 shows just the high $`B_2`$ tail of this distribution.

The probability distribution function, again a sum of delta functions, of $`h`$ the fraction of hydrophobic bits, is plotted in Fig. 7. As with Fig. 6, $`n_B=18`$ and $`ϵ=2`$. The distribution is peaked at $`h`$ a little above $`0.25`$: the mean value $`h=0.27`$ and the standard deviation is $`0.029`$. As $`h`$ increases towards $`0.5`$ then there are many more possible proteins but a rapidly increasing fraction of these are not soluble as a single phase according to our criterion. Thus there is a trade off between the number of possible proteins and the fraction that are soluble. This trade-off results in most proteins having between 20% and 35% of their bits hydrophobic. This of course depends on $`n_B`$ and $`ϵ`$. Increasing either one decreases $`h`$ but the picture remains qualitatively the same.

Let us return to our results for adenylate kinases, Fig. 3. Although it should be borne in mind that many of adenylate kinases’ hydrophobic amino acids will be in the centre of the protein, not at its surface, we can still quantify the scatter in $`h`$ for the kinases, and compare it to the scatter in $`h`$ for the model proteins. But of course any comparison will be purely qualitative. The adenylate kinases have around 200 amino acids in total, of which about 100 are classified as hydrophobic. We can try to model the distribution functions for $`h`$, for both adenylate kinases and our model proteins, with

$$h=n^1\underset{i=1}{\overset{n}{}}\zeta _i,$$
(5)
where for an adenylate kinase the sum is over all its amino acids in a protein, $`n=M`$, and for a model protein the sum is over the $`n=6n_B`$ bits. The $`\zeta _i`$ are independent random variables which are 1 with probability $`h`$ and zero otherwise. For adenylate kinases, see Fig. 3, we find that the standard deviation of $`h`$ is $`0.040`$, and Eq. (5) gives a standard deviation of $`0.035`$, only a little lower. To obtain the value of $`0.040`$ we took the sum over 206 terms; 206 is the mean length of the adenylate kinases in Fig. 3. Taking all the proteins to be the same length will decrease the spread slightly. Note that we can predict the distribution of the proteins’ hydrophobicity reasonably accurately using only the central limit theorem.

For our model proteins the standard deviation of $`h`$ is $`0.029`$, while Eq. (5) predicts $`0.043`$, which is rather larger but still comparable. Also, of course the shape of the distribution in Fig. 7 is quite close to Gaussian. Thus, the results for our model proteins are similar to those for real proteins, but as both are within a factor of $`1.5`$ of a simple prediction based on assuming the hydrophobic amino acids/bits are randomly distributed, it is hard to draw definite conclusions from this. The distribution of net charges can also be modeled assuming that the charged amino acids are distributed at random sear03 ; unpub .

## V Conclusion

We started with the idea that globular proteins needed to be soluble to function, and that their interactions depended on their surfaces which in turn were sensitive to which types of amino acids were at the surfaces of proteins. Then we defined a very simple model of a protein, whose surface-mediated-interactions depended on the values of strings of bits. A mutation in a protein such as an adenylate kinase which substituted a hydrophobic amino acid at the surface for a hydrophilic one could then be modeled by flipping one of these bits. Within our model, and with the constraint that a solution of the model protein is stable; the second virial coefficient is rather variable, its probability distribution function is plotted in Fig. 6. The criterion for the solution to be stable as a single phase is taken to be that the reduced second virial coefficient $`B_2/B_{2hc}1`$, which is enough for almost all fluids to be above their critical point. The condition that the protein solution be stable as a single phase is clearly a necessary condition, although in fact the second virial coefficient may be be more tightly constrained than this. Although the model used is simple, this variability does give credence to the idea that the variation in the fraction of hydrophobic amino acids in enzymes like adenylate kinase, see Fig. 7, gives rise to variability in the protein-protein interactions of these enzymes. In other words, that the second virial coefficients of E. coli’s and V. cholerae’s adenylate kinases may be significantly different, even though there is no obvious functional reason why their physical properties should differ. Unfortunately, virial coefficient measurements have not been performed for families of proteins. The variability is relevant to problems such as the purification and crystallisation of proteins. The separation of one protein from all the others in an extract from a cell which might contain thousands of proteins relies on differences in physical properties, charge, surface stickiness, etc., between proteins.

The probability distribution function of the second-virial coefficient, Fig. 6, is just the high $`B_2`$ tail of the distribution of all proteins. The remainder of the distribution function is cutoff by the requirement that $`B_2/B_{2hc}1`$. This full distribution function has a peak at an $`ϵ`$ dependent value of $`B_2`$; here well below $`1`$. Thus, without the cutoff at $`B_2/B_{2hc}=1`$, the distribution function is similar to the Gaussian distribution function found for the earlier model in which the patch-patch interactions were taken to be random variables sear02 . If we had kept with the previous model of describing with random variables the patch-patch interactions, and required that $`B_2/B_{2hc}1`$, then we would have obtained a distribution of second-virial coefficients similar to that in Fig. 6. In that sense a distribution like that in Fig. 6 is generic to any system where all model proteins except for those in a large $`B_2`$ tail are cutoff. However, within the earlier, simpler, model there is no clear way to look at either mutations and hence evolution, or to compare with sequence data for real proteins, as we did when we compared Figs. 3 and 7.

Finally, many simplifying assumptions have been made in order to arrive at our model system. It is therefore appropriate to comment on how this work can be extended to include more of the features of proteins inside cells. Both the model and our simple criterion for viability can be improved. The model is rather crude, and our sharp division between proteins deemed soluble and those deemed insoluble, could be softened. Then the fitness of a protein would decline over some range of values of the second virial coefficient. Also, we did not impose a maximum on the second virial coefficient. If it is important to limit the osmotic pressure, values of the second virial coefficient which are too positive may also be undesirable. However, in terms of understanding the behaviour of proteins in the complex crowded mixture of proteins that is the in vivo environment, perhaps the most important extensions of this work, is to multicomponent mixtures, and to include proteins which bind to each other. Inside cells thousands of different proteins are mixed together at a total protein concentration of around 20%, and many proteins are not monomeric but are part of complexes. The model studied here is flexible enough to both generate thousands of different proteins and to permit selective binding between proteins. Work on both is ongoing.

It is a pleasure to acknowledge that this work started with inspiring discussions with D. Frenkel. This work was supported by the Wellcome Trust (069242).

## Appendix: Computations

We are principally interested in the fraction of proteins that are soluble according to our criterion, and the distribution functions and means of various properties of soluble proteins. The fraction of proteins with $`B_2/B_{2hc}1`$ is determined by simply generating a very large number of proteins at random and finding the fraction that satisfy this requirement. The length of all runs are determined either by the requirement to obtain at least 2 significant figures or until longer runs produce almost identical plots. An exception is for $`n_B=18`$ and $`ϵ=2`$ where due to the smallness of $`f_v`$, it was only possible to obtain 1 significant figure of accuracy. The distribution functions, means etc., are obtained by starting with a soluble protein and generating a random walk in the space of soluble proteins. This is essentially no different from Metropolis Monte Carlo as applied to a system with a hard potential, e.g., a fluid of hard spheres, as our constraint $`B_2/B_{2hc}1`$, is a hard constraint. The averages are then obtained over these random walks.

The algorithm samples ‘protein space’ smith70 , with each sequence a unique point in this space and 2 sequences neighbours if 1 of them can be transformed into the other by a single mutation. This protein space is vast for real proteins and still very large for our model; for our model it contains $`2^{6n_B}`$ points. Note that all viable proteins are connected to all other viable proteins by an unbroken path of viable proteins and links between neighbouring viable proteins. This is easy to see if we consider that $`B_2`$ always either increases or stays the same if we flip a hydrophobic bit. Thus, starting from any viable protein we can flip each of its hydrophobic bits to hydrophilic bits, one at a time, until we reach the protein with all $`6n_B`$ bits hydrophilic. Each intermediate in this path must satisfy our solubility criterion as it is obtained from a protein which satisfies this criterion by flipping 1 or more hydrophobic bits. Thus we have proved that all viable proteins are connected to the protein with all hydrophilic bits, and so trivially all viable proteins are part of a connected network. This immediately implies that we can go from any one viable protein to any other via our Monte Carlo moves.


text:
# Comment on “Off-diagonal Long-range Order in Bose Liquids: Irrotational Flow and Quantization of Circulation”

## Abstract

In the context of an application to superfluidity, it is elaborated how to do quantum mechanics of a system with a rotational velocity. Especially, in both the laboratory frame and the non-inertial co-rotating frame, the canonical momentum, which corresponds to the quantum mechanical momentum operator, contains a part due to the rotational velocity.

A recent Letter su  claimed to have derived irrotational superflow or circulation quantization of a superfluid merely from the existence of off-diagonal long-range order (ODLRO). The derivation was based on an incorrect Hamiltonian.

In the beginning of that discussion, it is written: “Let us consider a bucket of Bose liquid composed of $`N`$ homogeneous interacting spinless particles, rotated with a constant angular velocity $`𝛀`$” su . This clearly indicates that $`𝐯_s=𝛀\times 𝐫`$ is the fixed rotation velocity of the bucket, not the superfluid velocity. $`𝐯_s`$ and $`𝛀`$ are both physical quantities, and one has no “gauge freedom” to add an additional term $`\theta (𝐫)`$ in $`𝐯_s`$, as appearing in Eq. (2) of su . With this understanding, the conclusion that $`𝛀`$ vanishes is, of course, absurd. In fact, for a superfluid in a rotating bucket, it is the Hamiltonian in the co-rotating frame that determines the equilibrium statistical mechanics landau ; pines . Instead of Eq. (1) of su , the relevant Hamiltonian is $`_j[\frac{(𝐩_jm𝛀\times 𝐫_j)^2}{2m}\frac{1}{2}m(𝛀\times 𝐫_j)^2]+\frac{1}{2}_{jl}V(r_{jl})`$, where $`𝐩_j=m(𝐯_j+𝛀\times 𝐫_j)`$ is the canonical momentum and is then replaced as $`i\mathrm{}_j`$, $`𝐯_j`$, $`𝐫_j`$ and $`𝐩_j`$ are relative to the reference frame co-rotating with $`𝛀`$ landau ; pines ; ll ; lq . This Hamiltonian obviously lacks the invariance under an $`𝛀\times 𝐫`$-dependent re-gauged translation underlying the argument employed later in su .

In their correspondence, the authors have claimed that they were not considering superfluid in a rotating bucket, instead, $`𝐯_s=𝛀\times 𝐫+\theta `$ was a “drift velocity” assumed to be a sum of an irrotational part $`\theta `$ and a rotational part $`𝛀\times 𝐫`$. Under this interpretation, even though the serious question about the physical meaning of this “drift velocity” is side-stepped, the Hamiltonian in su  is still incorrect. In fact, whenever a certain “drift velocity” contains a rotational part $`𝛀\times 𝐫`$, the canonical momentum relative to the non-inertial frame co-moving with this “drift velocity” is $`𝐩_j=m(𝐯_j+𝛀\times 𝐫_j)`$, and thus in the laboratory frame, the kinetic part of the Hamiltonian $`(𝐩_j+m\theta )^2/2m`$ ll , in contrast with Eq. (1) of su . It is very elementary knowledge that it is canonical momentum that corresponds to momentum operator in quantum mechanics ll ; lq . Moreover, in su , it is very mysterious why the assumed “drift velocity” suddenly became the superfluid velocity after the the rotational part was “shown” to vanish.

If the authors wish to adopt a reductio ad absurdam by showing that ODLRO is incompatible with the assumption that the superflow has a non-vanishing curl, i.e. the superfluid velocity is rotational, then the only valid interpretation is that $`𝐯_s`$ represents the superfluid velocity, first presumed to has a curl $`2𝛀`$; the purpose is then to show $`𝛀`$ has to vanish under ODLRO. If we use a subscript $`0`$ to represent the laboratory frame, then the assumption is $`\frac{1}{2}_0\times 𝐯_s=𝛀`$. With this interpretation, and if $`𝛀`$ is assumed to be position-independent, one can write $`𝐯_s=𝛀\times 𝐫_0+_0\theta (𝛀,𝐫_0)`$. In terms of the position vector $`𝐫`$ relative to the condensate, which moves with $`𝐯_s`$, the expression of $`𝐯_s`$ becomes very complicated. In general, Eq. (2) in su  is valid only if $`𝐫`$ and thus $`𝐩_j`$ in su  are relative to the rest frame. But this understanding cannot be consistent with the Hamiltonian in su . In fact, the Hamiltonian in the rest frame, in our notations, is simply $`_j\frac{𝐩_{j}^{}{}_{0}{}^{}{}_{}{}^{2}}{2m}+\frac{1}{2}_{jl}V(r_{jl})`$, where $`𝐩_{j}^{}{}_{0}{}^{}=i\mathrm{}_{j}^{}{}_{0}{}^{}`$, instead of Eq. (1) of su . If one subsequently obtains (1) by a “gauge” transformation, then $`m𝐯_s/\mathrm{}`$ is only the gradient of the phase factor, and $`𝐯_s`$ cannot be the superfluid velocity if it is presumed to be rotational as in (2) (cf. the preceding paragraph). On the other hand, it appears that $`𝐩_j`$ in su  are relative to the condensate and that the authors wish to give the Hamiltonian in the rest frame expressed in terms of the momenta relative to the condensate. With this interpretation, however, both Eq. (2) in su  and the expression $`𝛀=\frac{1}{2}\times 𝐯_s`$ below the equation are not valid generically under the presumption that the curl of the superfluid velocity is $`𝛀`$ in the laboratory frame. Therefore the Hamiltonian in su  is still incorrect.

It appears that the authors do not know that “irrotational superflow ” just means that the curl of the superfluid velocity vanishes, in the laboratory frame, and that the purpose they had set for their Letter was just to show that with ODLRO, the superfluid velocity is of the form $`\theta `$, i.e. its curl vanishes. Moreover, they ignore the well-known fundamental difference between rotation and translation landau ; pines ; ll ; lq . Note that in Sec. $`5.5`$ (pp. $`7576`$) of Ref. pines , it is a uniform (as italicized there) translation of a superfluid that is considered.

I thank Nigel Cooper for useful discussions.


text:
# Solitonic lattices in photorefractive crystals

## Abstract

Two-dimensional spatial solitonic lattices are generated and investigated experimentally and numerically in an SBN:Ce crystal. An enhanced stability of these lattices is achieved by exploiting the anisotropy of coherent soliton interaction, in particular the relative phase between soliton rows. Manipulation of individual soliton channels is achieved by use of supplementary control beams.

Wide ($``$1 mm) Gaussian beams launched in a photorefractive (PR) crystal in the self-focusing regime tend to break into spatially disordered arrays of filaments, owing to transverse modulational instabilities mamaev . However, ordered arrays of Gaussian beamlets ($``$10 $`\mu `$m), launched in conditions appropriate to the generation of spatial screening solitons opn02 , form much more stable solitonic lattices. Weakly interacting pixel-like arrangements of solitons that can individually be addressed are interesting for applications as self-adaptive waveguides petter ; chen .

Adaptive waveguides are of particular interest in all-optical information processing for their potential to generate large arrays, as well as for allowing many configurations with different interconnection possibilities. Spatial optical solitons are natural candidates for such applications, owing to their ability for self-adjustable waveguiding and versatile interaction capabilities, as demonstrated in light-induced Y and X couplers, beam splitters, directional couplers and waveguides. In addition to such few-beam configurations, the geometries with many solitons propagating in parallel—the so-called soliton pixels, arrays, or lattices—have been suggested for applications in information processing and image reconstruction mamyshev ; krol ; weilnau ; bram . Recently several groups demonstrated the formation of quadratic arrays of solitons in parametric amplifiers bram ; minardi  and PR media, with coherent weilnau ; petter2  and incoherent chen  beams.

In this communication we combine the properties of spatial PR solitons to form pixel-like lattices, to investigate experimentally and numerically the generation and interactions in large arrays of spatial solitons. We achieve improved stability of solitonic lattices by utilizing anisotropic interaction between solitons, in particular the phase-dependent interaction between solitonic rows. We manipulate individual or pairs of solitons using incoherent and phase-sensitive control beams.

Creation of solitonic lattices requires stable non-interacting propagation of arrays of self-focusing beams. A crucial feature in the parallel propagation of PR spatial solitons is their anisotropic mutual interaction krol2 . Because the refractive index modulation induced by a single soliton reaches beyond its effective waveguide, phase-dependent coherent as well as separation-dependent incoherent interactions, such as repulsion or attraction, may appear between the neighboring array elements. These interactions also affect the waveguiding characteristics of an individual solitonic channel. Therefore, the separation between solitons and their nearest-neighbor (NN) arrangement, determined by the form of the lattice, must carefully be chosen so as to minimize all forms of interaction. In our experiments, as well as numerical simulations, we determine the critical soliton distances for the given crystal thickness at which the interaction becomes noticeable. To achieve closer packing of solitons, while maintaining propagation without interaction, the phase relationship between different beams in the array is exploited.

The creation of lattices and controlling individual pixels is explored in the experimental setup of Fig. 1. Laser beam derived from a frequency-doubled Nd:YAG laser emitting at 532 nm illuminates a spatial light modulator (SLM), which imprints the image of a spot-array onto the beam. The spatial light modulator in turn is imaged onto the front face of a photorefractive SBN60:Ce crystal (5 x 5 x 20 mm), which is positioned so that the propagation direction is along the 20 mm axis. To exploit the dominant electro-optic component $`r_{33}200`$ pm/V of our SBN sample, the incident laser beam is linearly polarized parallel to the c-axis of the crystal, perpendicular to the propagation direction. Regular patterns of up to 25 x 25 spots, each with the diameter of about 20 $`\mu `$m and an intensity of 110 nW, are imaged onto the front face of the crystal. Applying an external DC electric field of about 1 kV/cm across the crystal and illuminating the crystal with uniform white light to induce an artificial dark conductivity, creates appropriate conditions for the formation of spatial screening solitons. To manipulate individual or pairs of solitons in the lattice, an additional control beam is derived from the same laser. It is our experience that large arrays of solitons can easily be formed in PR crystals. The larger the array the less significant the finite size (edge) effects. The number of pixels is mainly limited by the aperture of the PR crystal and the resolution of the inducing SLM. One can tailor the arrangement of beams according to specific needs.

Propagation without mutual interaction requires sufficient initial distance between beams. It also depends on other parameters influencing the soliton formation and propagation, such as the initial beam profiles and phases, and the strength of the nonlinearity, but less crucially. A distance smaller than the critical distance for coherent interaction between solitons leads to their attraction and eventual fusion. The attractive interaction is more pronounced along the y-transverse direction, perpendicular to the direction of the external field. For a given propagation distance this is reflected in the deformed appearance of the lattice. Deformations may also be caused by the crystal inhomogeneities that can not be controlled. Hence the determination of the critical distance is qualitative, based on the inspection of series of runs with decreasing inter-soliton distances. A typical series of runs is depicted in Fig. 2, where the distance between solitons is reduced by increasing the density of packing. An estimate of the critical distance is found to lie between 3 and 4 beam diameters for in-phase arrays in our experimental conditions. It is slightly larger in the $`y`$-direction than in the $`x`$-direction. This estimate is corroborated in numerical simulations.

Our numerical simulations are based on the paraxial approximation to the anisotropic propagation of coherent optical beams in saturable PR media,

$$2ikn_0_zA+_{}^2A=k^2n_0^4r_{33}(_x\varphi )A,$$
(1)
where $`k`$ is the wave number in vacuum, $`n_02.35`$ is the bulk refractive index, $`_{}^2`$ is the transverse Laplacian, and $`A`$ the slowly-varying envelope of the electric field of the beam. The paraxial equation is augmented with the equation for the electrostatic potential $`\varphi =𝐄`$ of the space-charge field $`𝐄`$ generated in the crystal zoz ; belic1

$$^2\varphi +\varphi \mathrm{ln}(1+I)=E_e_x\mathrm{ln}(1+I),$$
(2)
where $`E_e`$ is the external field. The light intensity $`I`$ is normalized to the background intensity. These equations are solved together, in the manner described in Ref. step . Behavior qualitatively similar to the experimental was found: Parallel propagation of coherent solitons with varying separation can be adjusted to be almost without interaction, (Fig. 3 (a)) whereas for separations below the critical length the fusion of soliton columns is observed. (Figs. 3 (b), (c)).

Propagation of incoherent lattices was found to be more stable than the propagation of coherent in-phase lattices for identical conditions, however our intention is to exploit the possibility of manipulating the distribution of phase across the lattice, to achieve closer packing while maintaining propagation without interaction. The comparison between incoherent, coherent and phase-engineered lattice is provided in Fig. 4. The NN separation is 40 $`\mu `$m, which is about half the critical distance from Fig. 2 and twice the FHWM of initial beams. It was found that the best results are obtained when the adjacent rows of solitons are out-of-phase. This is in agreement with theoretical expectations of the anisotropic interaction: Along rows the solitons are in-phase, and along columns they are out-of-phase, which makes the interaction between them universally repulsive and leads to enhanced stability. It is seen that the phase-engineered lattice in the right column of Fig. 4 is much more stable than the incoherent or coherent in-phase lattices in the left and central columns (note also different propagation distances). For smaller lattices, as in Fig. 4, the edge effects are also more pronounced.

To further test the idea of phase engineering for enhancing the stability, different lattices were formed and propagated (Fig. 5). It is seen that the in-phase lattices deform fast, and that after 10 mm of propagation, owing to the attractive forces along the y-direction, some rows are lost in the quadratic lattice, while the hexagonal and honeycomb lattices are strongly deformed. The honeycomb lattice deforms continually from the direct to the inverse (hexagonal) lattice and back (Figs. 5 (f) and (i)). However, out-of-phase lattices for identical initial conditions and at a longer propagation distance of 16 mm are only slightly deformed. It is also evident that the edge effects decrease with increasing the lattice size.

Use of solitonic arrays for applications in information technology requires means of manipulating individual waveguides, so as to combine different channels, separate them or induce energy exchange between the channels. For this purpose the well-known interaction effects of spatial PR solitons can be exploited opn02 ; weilnau . We utilize here a supplementary steering control beam derived from the Nd:YAG laser that can be focused anywhere on the front face of the crystal. Varying the phase of the control beam relative to the phase of array elements, a phase-sensitive coherent or incoherent interaction can be induced, that may lead to fusion or repulsion of different solitons in the array. Early experimental results on an incoherent control were presented in Refs. petter ; chen ; weilnau ; petter2 . Here we display the comparison with experiment, using an incoherent control beam acting on a pair of coherent pixels (Fig. 6). The data are as in the experiment of Refs. weilnau ; petter2 . Qualitative agreement is visible.

The action of a phase-sensitive control beam is depicted in Fig. 7. We chose to present the cases of 0 and $`\pm \pi /2`$ phase shifts of the control beam relative to the lattice. While for the zero phase shift the merging of the two channels is seen, in the case of $`\pm \pi /2`$ a strong energy exchange between the control beam and the solitons is observed. For $`\pi /2`$ phase shift the energy flows from the pixels to the control beam, so that at approximately 9 mm the control beam is maximally bright and the pixels are extinguished. For the phase shift of $`\pi /2`$ the flow of energy is from the control beam to the pixels. At approximately 10 mm the control beam is extinguished, and the two slightly displaced pixels are maximally bright. For the zero phase shift the three beams gradually merge, so that at approximately 15 mm only one beam is visible. Afterwards the cycles for the three values of the phase shift repeat in the inverse order, but less clearly visible, and some energy is transferred to the adjacent pixels (not shown).

In summary, we have investigated several aspects of the interaction of solitons arranged in lattices, for possible use in all-optical waveguides and interconnects. We have displayed experimentally the formation of large two-dimensional arrays of solitons that can be used to guide light at different wavelengths. Increased stability of solitonic lattices is achieved by the phase engineering of soliton rows. We have shown how to use a supplementary steering beam to fuse, extinguish, or enhance selected channels in the solitonic lattice. Our numerical simulations, based on the anisotropic model of PR solitons, are in qualitative agreement with experimental results.

AS and MB gratefully acknowledge financial support from the Alexander von Humboldt Foundation, for the stay and work at WWU Münster. DT acknowledges support from the Konrad-Adenauer-Stiftung. Parts of this work were supported by the Graduiertenkolleg ”Nichtlineare kontuierliche Systeme” of the Deutsche Forschungsgemeinschaft. Work at the Institute of Physics is supported by the Ministry of Science, Technologies, and Development of the Republic of Serbia, under grants OI 1475 and 1478. We thank Dr. Jürgen Petter, Institut für Angewandte Physik, TU Darmstadt, for help in the experiment on pixel control.


text:
# Effect of physical aging on the low-frequency vibrational density of states of a glassy polymer

## I INTRODUCTION

The effect of the physical aging on the macroscopic thermodynamical properties of glass are well known Struik (1978). On the other hand, the dependence of the structure at the nanometric scale on the physical aging is not yet clarified. It is expected that the possible effects on the nanostructure are very small and hardly observable by the diverse microscopy techniques. It was shown that the properties of vibrations, with energies in the meV spectral range, are dependent on the nanostructure because these vibrations are extended on nanometric lengths. For this reason, it is expected that the experimental study of the characteristics of low-energy vibrations will inform us about the glass nanostructure and its changes by physical aging.

It was earlier suggested that the excess of low-energy Raman scattering, called boson peak and the related excess of vibrational density of states (VDOS) (called also boson peak) is due to an inhomogeneous nanostructure: more cohesive nanodomains are separated by less cohesive zones Duval et al. (1990); Mermet et al. (1996a). From this model, it is expected that the contrast of bonding or of elastic constant between cohesive domains and softer zones weakens with aging, and that, in consequence, the vibrational density of states (VDOS) at low-energy decreases. From previous measurements, it is not clear that the aging induces a decrease of the VDOS at low-energy in the spectral range of the boson peak. Isakov et al. Isakov et al. (1993) showed, by inelastic neutron scattering, that the low-energy VDOS of the $`As_2S_3`$ glass increases after thermal quenching. This neutron measurement confirmed the increase of the low-temperature specific heat observed by Ahmad et al. Ahmad et al. (1986) after quenching in the same glass. On the other hand, no change of the VDOS was observed after quenching in the polybutadiene glass Zorn and Frick (1998).

Several experiments aimed to test the effect of thermal treatments on the VDOS have been carried out with the poly(methyl methacrylate) (PMMA) glass. Kanai et al. Kanaya et al. (1996) observed an increase of the inelastic incoherent neutron scattering (INS) intensity at low-energy after quenching from a temperature of $`453K`$ of an annealed sample. However, as it will be shown in the discussion, the conditions of this experiment do not allow to conclude that the rejuvenation (by thermal quenching) and the physical aging (by thermal annealing at a temperature lower than the glass transition temperature, Tg) have opposite effects on the low-energy VDOS. New experiments aiming to test the effect of aging on the low-energy VDOS of PMMA were carried out recently Etienne et al. (2002). It was found that the aging by annealing at $`363K`$, i.e. $`10K`$ below Tg, of a slightly cross-linked PMMA has no effect on the VDOS determined from INS at low temperature, $`T=30K`$. On the other hand, by methanol-assited aging at room temperature of the same PMMA Etienne et al. (2001) a significant decrease of the VDOS in the spectral range of the boson peak was observed Etienne et al. (2002). These observations were in agreement with the Raman ones. A decrease of the Raman boson peak was observed after methanol-assisted aging of the cross-linked PMMA at $`300K`$, and none after annealing at $`363K`$. However, the effect of the methanol-assisted aging was much stronger with non cross-linked PMMA Etienne et al. (2001). Probably, the cross-linking makes the nanostructure change much more difficult by annealing.

From the described situation, it is not yet clear that the low-energy VDOS of non cross-linked PMMA decreases by aging, i.e., by annealing at a temperature slightly lower than Tg. A decrease of the Raman boson peak at low-temperature was observed. But the low-energy Raman scattering is proportional to the VDOS and to the light-vibration coupling coefficient C(E), and it is possible that this is the decrease of this coefficient which is responsible for the decrease of the Raman boson peak and not the VDOS. To solve this problem new measurements of inelastic neutron scattering at low-temperature were carried out to determine directly the effect of the physical aging on the VDOS of PMMA. The comparison of the VDOS of an aged PMMA with the one that is rejuvenated by thermal quenching is presented in this paper. The new neutron experimental results with others obtained earlier are discussed to interpret the effect of aging on the glass nanostruture.

## II EXPERIMENTAL

The PMMA plates of clinical grade without any additives were purchased from Goodfellow SARL (reference ME303031). The average molecular weight is close to $`10^6`$ and the glass transition temperature $`T_g`$ is about $`390K`$, as determined by differential scanning calorimetry (DSC). PMMA disks with a diameter of 5 cm and a thickness of about 0.04 cm were cut for heat treatments and inelastic neutron scattering.

The inelastic neutron spectra were recorded on the time-of-flight instrument IN6 at ILL, Grenoble. The energy resolution is $`\delta E=80\mu eV`$. The momentum transfer range extends from $`Q=0.22\text{Å}^1`$ to $`Q=2.06\text{Å}^1`$. The spectra were taken in the temperature range $`2300K`$, using an helium cryofurnace. The temperature of $`30K`$ was chosen to obtain the neutron inelastic scattering by the low-energy harmonic vibrations much higher than the one by anharmonic or relaxational motions Mermet et al. (1996a). The scattering cross-sections were obtained after the usual standard calibrations by means of the vanadium runs and the removal of the empty-can contributions. The VDOS for harmonic modes were obtained by taking the average of the spectra given by the different detectors, that is the average over the range from $`Q=0.22\text{Å}^1`$ to $`Q=2.06\text{Å}^1`$. It was calculated through the use of an iterative procedure described elsewhere Fontana et al. (1990). The so-obtained VDOS were corrected by the Debye-Waller factor and the multiphonon contributions.

## III EXPERIMENTAL RESULTS

Two PMMA samples were compared. Sample-1 was heated for 30 minutes at $`413K`$, quenched at room temperature, annealed for 30 days at $`373K`$ and then 15 days at $`363K`$. In the view of the energy landscape, the annealing at two different temperatures was thought to be more efficient for aging: while the higher temperature allows the system to cross over high energy barriers, the lower temperature allows a stabilization into deeper wells. Sample-2 was heated for 30 minutes at $`413K`$ and quenched at room temperature on a copper plate. The cooling rate was estimated to be between 50 and 100 K/s. It was controlled, by the observation of the Raman lines, that these thermal treatments, especially the heating at $`413K`$, did not induce depolymerization. It was checked that the intensity of the line at $`1640cm^1`$, that corresponds to the vibration localized on the C=C bond, is zero for both samples. It means that the concentration of MMA monomers is very weak, and that there is no depolymerization Surovtsev et al. (1996) even after heating at $`413K`$ for 30 min. As a consequence, Sample-1 was physically aged and Sample-2 was rejuvenated, its thermal history being erased Struik (1978).

The effect of aging was clearly observed by DSC. A strong endothermal peak appeared for Sample-1, and not for Sample-2, as one can see in Fig. 1. In addition, it was determined, from the DSC curves, that the fictive temperatures were respectively $`T_f=366K`$ and $`T_f=390K`$ for Sample-1 and Sample-2, in agreement with the thermal treatments.

The incoherent scattering functions of both samples, $`S(E)`$, taken at $`T=30K`$, are compared in Fig. 2. They were determined by taking the average of the spectra given by the different detectors over the range from $`Q=0.22\text{Å}^1`$ to $`Q=2.06\text{Å}^1`$. They are corrected for the elastic peak, that was obtained by measuring the inelastic scattering from the samples at $`T=2K`$. As their shapes are identical from the energy equal to 4 meV, they were normalized by coincidence of the curves from the energy of 4 meV for a better comparison. Considering the experimental error bars due to neutron counting, which are shown on the curves at 1.25 meV, it is clear that $`S(E)`$ is more intense in the $`0.62meV`$ spectral range for Sample-2 than for Sample-1.

The VDOS $`g(E)`$ were deduced from $`S(E)`$ Fontana et al. (1990). These VDOS divided by the square of energy, $`g(E)/E^2`$, are plotted in Figure-3. As expected from $`S(E)`$ (Figure-2), the boson peak at 1.75 meV is more intense for the rejuvenated Sample-2 than for the aged Sample-1, even if the difference is not huge. Furthermore, it is evident, from Figure-3, that the observed difference originates from the VDOS for harmonic vibration modes and not for anharmonic ones, because it does not decay monotonously from E=0 to higher energy as expected for anharmonic motions: it is observed (inset of Figure-3) that this difference presents a maximum around E = 1.05 meV. It was detected that it is not the case for a sample slightly depolymerized by heating at 413 K for a much longer time and containing some monomers: even at $`T=30K`$ neutron scattering by anharmonic motions due to the presence of monomers was observed, in agreement with previous Raman measurements Surovtsev (1996). The VDOS were also determined at $`T=300K`$. However, at this temperature, the neutron scattering by anharmonic motions is relatively high for both samples and masks the scattering by harmonic vibrations at low-energy.

## IV DISCUSSION

In this described experiment much care was taken, on the one hand, to make sure that the heat treatments induce only physical aging and not a chemical one like depolymerization, and, on the other hand, to obtain the effect of aging on the VDOS corresponding to harmonic vibration modes. In a previous study of the effect of annealing on the low-energy excitations in PMMA, the conditions of the experiment were different Kanaya et al. (1996). Indeed, these authors of this study compared a sample annealed 7 days at $`T=363K`$ to a sample quenched from $`T=453K`$. It was shown that at this temperature depolymerization occurs and induces an increase of the monomer concentration Surovtsev et al. (1996). The reversible effect was observed by subsequent annealing at a temperature lower than Tg Surovtsev et al. (1996). As a consequence, by such a thermal treatment Kanaya et al. (1996), this is a partially chemical aging which is obtained instead of a pure physical aging. Furthermore, in this experiment Kanaya et al. (1996) the inelastic scattering was measured at $`T=73K`$. At this temperature the scattering by anharmonic or relaxational motions is not negligible Mermet et al. (1996a). In these conditions and after observation of the published spectra, it is likely that the difference in the inelastic neutron scattering of annealed and quenched PMMA observed by Kanaya et al. is mainly due to anharmonic motions.

In another glassforming polymer, polybutadiene, no effect of annealing was observed on the VDOS by measuring the inelastic neutron scattering Zorn and Frick (1998). It could be suggested that the decrease of the low-energy VDOS by aging is not general for glass-forming polymers. However, the used thermal treatments can explain the non-observation of difference in the VDOS between annealed and quenched polybutadiene. The aged sample was heated for 36 h at $`T=169K`$ (11 K below Tg), and the quenched sample was heated for 1 h at 270 K and cooled at 100 K with a rate of 2 K/min. The cooling rate used for this quenching was possibly too slow to avoid an undesidered effect of aging during the cooling from 270 K to 100 K.

In a very recent study of a crystallisable silicate material by inelastic neutron scattering, it was observed that the intensity of the boson peak of a hyperquenched ($`10^6K/s`$) sample decreased strongly after annealing at a temperature below Tg Angell et al. (2003). The neutron measurements were performed at room temperature, and probably the measured VDOS is partially accounted for anharmonic vibrational modes. A similar behavior was earlier observed with a metallic glass Suck (1989). These results were discussed in the frame of the potential energy landscape model Angell et al. (2003).

A simple explanation of the decrease of the VDOS by aging, in our experiment, could be the increase of sound velocity concomitant with the increase of density, and in consequence a decrease of the Debye VDOS contribution Isakov et al. (1993). However, the difference of $`g(E)/E^2`$ between aged and rejuvenated PMMA is expected to be constant as a function of E. It is not the case, this difference has a maximum around 1.05 meV (inset of Figure-3). Instead of having an effect on the Debye VDOS, the aging would change the excess of VDOS responsible for the boson peak. The interpretation of the low-frequency modes in the boson peak is still a matter of debate. Several theoretical models are based on the system dynamics without consideration of the structure, among which the very recent landscape model Grigera et al. (2003). Here, we will consider the nanostructural point of view, which is not in contradiction with the dynamical approach and may be relevant to discuss the aging effect. The boson peak was previously Duval et al. (1990); Mermet et al. (1996a) interpreted by an inhomogeneous glass elasticity or cohesion at the nanometric scale: more cohesive domains would be separated by softer zones. This view is supported by a recent simulation of the elasticity of small amorphous systems Wittmer et al. (2002). Previous experimental results are in agreement with quasi-periodic arrangements of the cohesive domains on short distances Mermet et al. (1998); Duval and Mermet (1998), so that the boson peak would correspond to a Van Hove singularity associated to the quasi-periodic arrangements Duval et al. (2002). Such interpretation can be compared to the model of Taraskin et al. Taraskin et al. (2001), in which the boson peak is related to the lowest Van Hove singularity of the reference crystalline system pushed down by the effect of disorder-induced level repulsion. In the frame of the inhomogeneous glass nanoelasticity model, the decrease of the low-energy VDOS by aging can originate from a decrease of the contrast in the elasticity at the nanometric scale. This would be in agreement with the observations of Isakov et al. about the effect of aging in the $`As_2S_3`$ glass Isakov et al. (1993). These authors determined that the relative effect of aging on the low-energy VDOS was the same as the one on the light-vibration coupling coefficient $`C(E)`$ in Raman scattering. This coefficient $`C(E)`$ increases with the amplitude of the dielectric susceptibility fluctuations Jäckle (1981). The dielectric susceptibility is directly related to the elastic constants by the elasto-optic tensor. Therefore, the dielectric susceptibility contrast at the nanoscale goes like the elastic constant one, and one may conclude that the decrease of $`g(E)`$ and $`C(E)`$ have a common origin: the smoothing of the elastic constant contrast at the nanometric scale. The same conclusion can be proposed for PMMA. Unfortunately, it was not possible to perform a quantitative comparison between $`g(E)`$ and $`C(E)`$ in the case of PMMA, because the effect of aging on these quantities are relatively much weaker for this glass than for $`As_2S_3`$ Ahmad et al. (1986); Isakov et al. (1993).

The results obtained by the inelastic neutron scattering study presented in this paper must be compared to other previous ones. Same measurements were carried out with PMMA slightly cross-linked by trimethyl-1,1,1-propane trimethacrylate (TRIM) Etienne et al. (2002). With this cross-linked PMMA, no effect of aging was observed on the low-energy Raman scattering by harmonic vibration modes. It is believed that it is due the presence in the zones between domains of cross-links which reinforce the linking of a domain with its next-neighbors and, then, precludes the decrease of the contrast between soft zones and cohesive domains with aging. Methanol-assisted aging of PMMA at room temperature was also studied Etienne et al. (2001, 2002). PMMA plates was soaked for more than one month in methanol at room temperature and then dried. By mechanical relaxation measurement, it was shown that the effect of this treatment corresponds to an aging at room temperature accelerated by the presence of methanol. This can be explained by a preferential location of the methanol in the softest zones, which makes easier the molecular segmental motions.

It may be noticed that the tendency towards a more homogeneous structure by aging was deduced from elastic light scattering by Takahara et al. in the case of PMMA Takahara et al. (1999). These authors determined the correlation length and the amplitude of density fluctuations at the scale of 200 nm. They found that the values of these two fluctuation characteristics decreased after an aging at $`T=353K`$ for less than 10 hours. These conclusions are not exactly comparable to the ones suggested in this paper. In the present work, the explored length scale is 10 times smaller, and deals with the fluctuations of elastic constant, and not with the density ones.

The effect of aging or rejuvenation by thermal treatments on the low-energy vibrations, is now compared to that of a plastic deformation Mermet et al. (1996b). It is well-known that a plastic deformation rejuvenates a glass Struik (1978); Hasan and Boyce (1993); Hasan et al. (1993); Hodge (1995); Utz et al. (2000). The effect of plastic deformation by shearing on the low-energy vibrations was investigated by Raman scattering in the case of PMMA and polycarbonate (BPA-PC) glasses Mermet et al. (1996b). An increase of the Raman scattering intensity in the region of the boson peak was observed. This increase of Raman intensity can be attributed, at least partially, to an increase of VDOS. Furthermore, as a confirmation, no increase of low-energy Raman scattering was observed after plastic deformation of the cross-linked PMMA Mermet et al. (1996b), like after thermal quenching Etienne et al. (2002). These results confirm the effect of aging (and rejuvenation) on the low-energy VDOS and on the glass nanostructure, consistently with the proposed interpretation.

It is noted that the experimental results and their interpretation are in agreement with, on the one hand, the view of a vitrifying liquid, that was deduced by Stillinger Stillinger (1988) from the potential energy landscape, ”as a patchwork of relatively strongly bonded molecular domains separated by irregular walls of weakened bonds”; and, on the other hand, with the recent numerical simulation of aging and rejuvenation, which shows that the aging induced a pronounced decrease in the system pressure, and that local compression and expansion of rearranging clusters of atoms appear as the shear deformation proceeds and, therefore, erases the aging Utz et al. (2000). This description is not far from the one involving the arrangement of more cohesive domains separated by less cohesive zones as proposed earlier Duval et al. (1990); Mermet et al. (1996a, 1998); Duval and Mermet (1998), the contrast between cohesive domains and softer zones decreasing with aging, and inversely increasing with rejuvenation, as suggested in this discussion.

Finally, it is remarked that the observed decrease of the VDOS with aging confirms qualitatively the very recent predictions deduced from the energy landscape model Grigera et al. (2003).

## V CONCLUSION

The inelastic neutron scattering experiment described in this paper shows clearly a lower harmonic vibration density of states for aged poly(methyl methacrylate) than for thermally quenched (rejuvenated) one. Much care was taken to avoid chemical effects by thermal treatment. The neutron measurements were carried out at low-temperature to obtain the scattering by harmonic vibrations. The effect of thermal quenching is compared to the one of a plastic shearing observed earlier. The rejuvenations obtained by both thermal and mechanical treatments are characterized by an increase of the low-energy vibrational density of states. These effects on the vibrational density of states are in agreement with the model of an inhomogeneous cohesion in the polymeric glasses at the nanometric scale, and with the energy landscape deductions. It would be interesting to establish a relation between the glass nanostructure or the nanoelasticity and the profile of the energy landscape.

###### Acknowledgements.

The authors are grateful to J. Ollivier from ILL for his assistance in the neutron measurements, and to A. Mermet for discussion. They thank warmly C. A. Angell to send them a copy of his paper Angell et al. (2003) before publication and for his comments.


text:
# 1 Introduction

## 1 Introduction

If, according to the assumptions of statistical physics, the second law is regarded as a “fact” rather than a dynamical law, it could conceivably not hold at all, hold only occasionally, or even apply in varying directions of time. Larry Schulman demonstrated very nicely and convincingly in his oral contribution to this conference how this may happen in principle (see also ). The major remaining question is whether his examples can be regarded as realistic in our universe.

In particular, we may understand from Schulman’s examples how a certain time arrow depends on “improbable” (low entropy) initial or final conditions – regardless of the direction in which we perform our calculation. The latter (apparently trivial) remark may be in place, since many derivations of the second law tacitly assume in a crucial way that the calculation is used to predict. That is, it is assumed to follow a “physical” direction of time (from an initially given present towards an unknown future). However, precisely this physical arrow, or the fact that only the past can be remembered and appears “fixed”, is a major explanandum.

While, for a given dynamical theory, we know in general precisely what freedom of choice remains for initial or final conditions, mixed ones (such as two-time boundary conditions) are subject to dynamical consistency requirements – similar to an eigenvalue problem with given eigenvalue. This problem remains relevant even for incomplete (for example, macroscopic) initial and final conditions. It is usually difficult to construct an individual solution that is in accord with both of them. In Schulman’s examples, individual solutions were mostly found by “trial and error”, that is, by exploiting a sufficient number of solutions with given initial conditions and selecting those which happen to fulfill the final ones (or vice versa). However, in a realistic situation it would be absolutely hopeless in practice ever to end up with the required low entropy because of the exponential growth of probability with entropy. Only an exponentially small fraction of all solutions satisfies one or the other low entropy boundary conditions. In the case of complete mixing, it is the square of this very small number that measures the fraction of solutions with two-time boundary conditions.

Being able to find solutions by trial and error thus demonstrates already the unrealistic case. This difficulty in finding solutions does not present any problem for their existence on a classical continuum of states if mixing is sufficiently complete: any set of solutions with finite measure can be further partitioned at will, since entropy has no lower bound in this classical situation. This conclusion is changed in quantum theory, which would in a classical picture require the existence of elementary phase space cells of Planck size $`h^{3N}`$. The product of initial and final probabilities characterizing the required low-entropy may then represent a phase space volume smaller than a Planck cell – thus indicating the absence of any solution.

## 2 Retarded and Advanced Fields

The consistency problem in a classical setting (though without mixing) is discussed in Schulman’s “wet carpet” example, intended to prove the compatibility of two interacting systems with different, retarded and advanced, electrodynamics . It is similar to an example studied by Wheeler and Feynman , where a charged particle, bound to pass an open trap door, is assumed to shut it before the particle’s arrival by means of advanced fields if it ever enters the space behind that door. In both examples, there is but a very narrow band of consistent solutions, in Schulman’s example represented by a partly opened window. This narrow band may exist for systems which are far from thermodynamical mixing (cf. the following sections), and they may be consistent only if the model is considered in isolation. In reality, macroscopic objects always interact with their surroundings. In a causal world, this would produce “consistent documents” (not only usable ones) in the thermodynamical future. In this way, information may classically spread without limit, thus leading to inconsistencies with an opposite arrow of other systems. In quantum description, classical concepts even require the presence of irreversible decoherence, while microscopic systems would remain in quantum superpositions of all their histories (see Sect. 6).

Philosophers are using the term “overdetermination of the past” to characterize this aspect of causality . (Note that the conventional additive physical entropy neglects such nonlocal correlations, which would describe the consistency of documents, for being dynamically irrelevant in the future .) In a deterministic world, one would thus have to change all future effects in a consistent way in order to change the past. For example, classical light would even preserve its in principle usable memory forever in a transparent universe. In classical electrodynamics, there is but one real Maxwell field, while the retarded and advanced fields of certain sources are merely auxiliary theoretical concepts. The same real field can be viewed as a sum of incoming and retarded, or of outgoing and advanced fields, for example (see Chap. 2 of ). Retarded and advanced fields (of different sources) thus do not add. Observing retarded fields (as our sensorium and other registration devices evidently do) means that incoming fields related to unspecified past sources (“noise”) are negligible – incompatible with the presence of distinctive advanced radiation.

Problems similar to those with opposite arrows occur with closed time-like curves (CTCs), which are known to exist mathematically in certain solutions of Einstein’s field equations of general relativity. This existence means that local initial and final conditions for the geometry, defined with respect to these closed time-like curves, are identical and thus dynamically consistent. However, CTCs are incompatible with an arrow of time, such as an electrodynamic or thermodynamic one. Those clever science fiction stories about time travel, which are constructed to circumvent paradoxes, and thus seem to allow CTCs for human adventurers, simply neglect all irreversible effects which must arise and would destroy dynamical consistency. Since geometry and matter are dynamically coupled, boundary conditions which lead to an arrow of time must also protect chronology (whatever the precise dynamical model). Wet carpet stories belong to the same category as science fiction stories: they do not resolve the unmentioned paradoxes that would necessarily arise from opposite arrows.

## 3 Cat Maps

Borel demonstrated long ago  that microscopic states of classically described gases are dynamically strongly coupled even over astronomical distances. This is a time-symmetric consequence of their extremely efficient chaotic behavior, caused by deterministic molecular collisions. Of course, this does not mean that macroscopic properties are similarly sensitive to small perturbations, although their fluctuations (such as Brownian motion) must be affected.

Macroscopic properties characterize the microscopic state of a physical system incompletely, for example by representing a coarse graining in phase space (or, more generally, a Zwanzig projection – Chap. 3 of ). The deterministic dynamics of initially given coarse grains is often described by measure-preserving dynamical maps. In contrast to deformations of extended individual objects in space (such as Gibbs’ ink drop), and even in contrast to the N discrete points in single-particle phase space which represent a molecular gas, Kac’s symbolic “cats” (areas in phase space)  represent ensembles, or sets, of possible physical states of a given system. Therefore, Schulman’s entropy  as a function of deformed cats (“cat maps”) is an ensemble (or average) entropy – not the entropy of an individual physical state. The entropy of this ensemble is defined to depend on its distribution in phase space, obtained after coarse graining with respect to given and fixed grains as a macroscopic reference system, while the entropy of an individual state (point in phase space) would be given solely by the size of the specific grain that happens to contain it at a certain time.

This is essential (and sufficient) for Schulman’s argument that the intersection of two sets representing specific initial or final conditions is not empty if mixing is complete. In our universe, however, some variables participate in very strong mixing, while others (“robust” ones, such as electromagnetic waves or atomic nuclei) may remain stable for very long times. They are the ones that may store usable information.

Since cat maps describe sets of states for rather simple dynamical systems, their dynamics is far less sensitive to weak interactions than that of individual Borel type systems. For this reason, two systems described by cat maps with opposite arrows of time may even be consistent for mild interactions . However, these cat maps do not form a realistic model appropriate to discuss thermodynamical arrows in our universe.

## 4 (Anti-)Causality

In order to define causality without presuming a direction of time, one has to refer to the internal structure of the evolving dynamical states. The above-mentioned over-determination of the past (in other words, the existence and consistency of multiple documents) is a typical example. Another one is given by the concentric waves emitted from a local source. In our world, both are empirically (not logically) related to a time direction.

While one may expect that all such internal structures can be shown to evolve in time from appropriate initial conditions, they are too complex to be investigated in terms of Schulman’s simple models. For example, retarded (concentrically outgoing) waves exist in the presence of sources precisely when incoming fields are negligible. This can be the case “because” of an initial condition for the fields, or because of the presence of thermodynamic absorbers .

Instead of these specific structures in the states, Schulman studied the “effect” (in both directions of time) of “perturbations” defined by small terms changing the Hamiltonian at a certain time only . This effect is not easily defined in a time-symmetric way, since an “unperturbed solution” defined on one temporal side of the perturbation would be exclusively changed on the other one (no matter which is the future or past). If the unperturbed solution obeyed a two-time boundary condition, the perturbed one would in general violate it on this “other” temporal side. In contrast to the above mentioned internal structures, our conventional concept of perturbations is based on the time direction used in the definition of external operations.

Therefore, in a first step, Schulman considered sets of solutions again. The set of all solutions obeying the “left” boundary condition (in time) remains unchanged on the left of the perturbation, while the opposite statement is true on the right. However, individual solutions found in the intersection of these two sets (consisting of those ones which fulfill both boundary conditions) in the case of a perturbation are generically different from the unperturbed ones on both sides of the perturbation. Now, if mixing is essentially complete on the right of the perturbation (that is, for a sufficiently distant right boundary), the right boundary condition does not affect the solutions which form the intersection (considered as a set) on the left. This means that mean values of macroscopic variables in the set of all solutions that are compatible with both boundary conditions may only differ on the right (a consequence regarded as causality by Schulman) . This is true, in particular, for the mean entropy (if the latter is defined as a function of macroscopic, that is coarse-grained, variables). Individual solutions can not be compared in this way, since there is no individual relation between them. In the case of complete mixing on the right, there is even a small but non-empty subset of solutions of the original two-time boundary value problem which keep obeying the right boundary condition without being changed on the left. However, using them for the argument would mean that only very specific solutions can be perturbed in this specific sense.

In a second approach, Schulman studies the “effect” of macroscopic perturbations on individual solutions of an integrable system. This system is defined as consisting of a finite number of independent oscillators with different frequencies. Although solutions which fulfill both boundary conditions can be found with and without an appropriate perturbation, they are again not individually related. Therefore, the causal interpretation of the perturbation remains obscure. (For closed deterministic systems, any perturbation would itself have to be determined from microscopic boundary conditions, and the consistency problem becomes even more restrictive than for just two boundary conditions.)

Nonetheless, I was pleased to discover that Schulman’s model is formally identical with a model of particles freely moving on a periodic interval (a “ring”) that I had used in an appendix of  for much larger numbers of constituents than used by him (such that finding two-time boundary solutions by trial and error would be hopeless). Particle positions on the ring have merely to be re-interpreted as oscillator amplitudes in order to arrive at Schulman’s picture. I used this opportunity to search by trial and error among analytically constructed two-time boundary solutions for those ones which happen to possess slightly lower entropy than the mean at some given “perturbation time” $`t_0`$ (see Fig. 1). Unfortunately, the results do not confirm Schulman’s claim that these solutions are “affected” by the perturbation only in the direction away from the relevant low entropy boundary (that is, towards the “physical future”) . Evidently, this concept of causality, defined by means of perturbations, is insufficient. The very concept of a “perturbation” seems to be ill-defined for two-time boundary conditions.

As another example, I calculated the effect on the solution in both directions of time that results from a microscopic perturbation of the state (in this case simply defined by an interchange of velocities between particles at some time $`t_0`$). Both boundary conditions are then violated by the new solution arising from this perturbed state, used as an “initial” condition. The results (shown in Fig. 2) are now most dramatic towards the former “past”, demonstrating the relevance of fine-grained information for correctly calculating “backwards in time”. Deviations from the original two-time boundary solution close to $`t_0`$ also on the right are due to the fact that the coarse graining assumed in this model does not define a very good master equation (as discussed in ).

## 5 Cosmology and Gravitation

It appears evident from the above discussion that opposite arrows of time would in general require almost complete thermalization between initial and final conditions, which is hard to accomplish even in a cosmological setting. In its present stage, this universe is very far from equilibrium. A reversal of the thermodynamical arrow of time together with that of cosmic expansion, as suggested by Gold , would therefore require a total life time of the universe vastly larger than its present age. Weakly interacting fields may never thermodynamically “mix” with the rest of matter .

In particular, Davies and Twamley demonstrated  that an expanding and recontracting universe would remain essentially transparent to electromagnetic waves between the two radiation eras. This means that advanced radiation resulting from all stars which will exist during the recontraction of our universe would be present now, apparently unrelated to any individual future sources because of their distance, but red- or blue-shifted, depending on the size of the universe at the time of (time-reversed) emission. According to Craig , this radiation would show up as a non-Planckian high-frequency tail of the cosmic background radiation resulting from the past radiation era (where it would be absorbed in time-reverse description). This leads to the consistency problems described in Sect. 2.

While neutrinos from the future would presumably remain unobserved, gravity, despite its weakness, dominates the entropy capacity of this world, and leads to consequences which are the most difficult ones to reverse. Black holes are expected to harbor event horizons which would not be able ever to disappear in classical relativity, while in quantum field theory they are predicted to disappear into Hawking radiation in the distant future in an irreversible manner. However, contraction of gravitating objects, including the formation of black holes, requires that higher multipoles are radiated away. This radiation arrow is the basis of the “no hair theorem”, which would characterize the asymptotic final states of black holes in an asymptotically flat and time-directed universe. Because of the diverging time dilation close to a horizon, any coherent advanced radiation (with the future black hole as its retarded cause) would be able to arrive in time to prevent the formation of an horizon . This solution of the “information loss paradox” may save a deterministic universe (without leading to inhomogeneous singularities).

While the required initial and final conditions are not obviously consistent in this classical scenario, this problem is relaxed in quantum cosmology.

## 6 Quantum Aspects

Realistic models of physical systems require quantum theory to be taken into account. Since quantum entropy is calculated from the density matrix (that may result from a wave function by means of generalized “coarse graining”), its time dependence has in principle to include a collapse of the wave function during measurements or other “measurement-like” situations (such as fluctuations or phase transitions). If the collapse represents a fundamental irreversible process, it defines an arrow of time that is never reversed. Only a universal Schrödinger equation (leading to an Everett interpretation) could be time (or CPT) symmetric. A reversal of the time arrow would then require decoherence to be replaced by recoherence: advanced Everett branches must combine with our world branch in order to produce local coherence. Although being far more complex than a classical model (since relying on those infamous “many worlds”) this would still allow us to conceive of a two-time boundary condition for a global wave function (see Sect. 4.6 of ). Note that Boltzmann’s statistical correlations (defined only for ensembles) now become quantum correlations (or entanglement, defined for individual quantum states). For example, re-expanding black holes, mentioned in the previous section, would in an essential way require (and possibly be facilitated by) recoherence.

In quantum gravity (or any other “reparametrization-invariant” theory), the Schrödinger equation is reduced to the Wheeler-DeWitt equation, $`H\mathrm{\Psi }=0`$, which does not explicitly depend on time at all. However, because of its hyperbolic form, this equation defines an “intrinsic initial value problem” with respect to the expansion parameter $`a`$. In a classical (time-dependent) picture, the initial and final states would have to be identified in order to define one boundary condition, while the formal final condition (with respect to $`a`$) for recontracting universes is reduced to the usual normalizability of the wave function for $`a\mathrm{}`$. Big bang and big crunch (distinguished by means of a WKB time, for example) could not even conceivably be different as (complete) quantum states (Chap. 6 of ), while forever expanding universes might be said to define an arrow of time that never changes direction during a WKB history. All arrows thus seem to be strongly entangled.