text:
# Viscosity in the escape-rate formalism

## I Introduction

The link between the irreversible phenomena governed by macroscopic equations such as the Navier-Stokes equations and the microscopic reversible dynamics of the atoms and molecules is a fundamental problem. In this context, it has been shown that typical many-body systems of interacting particles present a chaotic dynamics livi ; posch-hoover ; vanbei-dorfman ; gasp-vanbei . This microscopic chaos develops a sensitivity to initial conditions over a time scale of the order of the intercollisional time of the atoms and molecules. This sensitivity to initial conditions is characterized by positive Lyapunov exponents which are the rates of exponential separation between some reference and perturbed trajectories of the system. The sensitivity to initial conditions results into a huge dynamical randomness characterized by a positive Kolmogorov-Sinai (KS) entropy per unit time given by the sum of positive Lyapunov exponents if the system is at equilibrium:

$$h_{\mathrm{KS}}=\underset{\lambda _i>0}{}\lambda _i,$$
(1)
an identity known as Pesin’s equality eckmann-ruelle ; ott . This dynamical chaos provides an efficient mechanism of randomization of the different observable quantities such as the microscopic currents associated with the transport properties.

For a Hamiltonian-type microscopic dynamics, a connection can be established between dynamical chaos and the transport properties thanks to the escape-rate formalism gasp-nicolis ; dorf-gasp ; gasp-dorf ; gasp-book ; dorf-book . In this formalism, the gap between the kinetic time scale of chaotic properties and the hydrodynamic time scale of transport properties is bridged by linking the transport coefficients to differences between chaotic properties. In this formalism, the following formula has been derived for viscosity dorf-gasp

$$\eta =\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2\left(\underset{\lambda _i>0}{}\lambda _ih_{\mathrm{KS}}\right)__\chi ,$$
(2)
where the difference between the sum of positive Lyapunov exponents and the KS entropy is nonvanishing because the system is here under nonequilibrium conditions. These nonequilibrium conditions select the trajectories of the many-body system which do not escape out of a phase-space region specific to the transport property of interest. This phase-space region is defined by requiring that the Helfand moment associated with the transport property remains bounded in an interval of extension $`\chi `$. In the limit $`\chi \mathrm{}`$, the nonequilibrium condition is progressively relaxed and the sum of positive Lyapunov exponents as well as the KS entropy approach their equilibrium value satisfying Pesin’s equality (1). Under nonequilibrium conditions, Pesin’s equality is not satisfied and the difference gives the rate of escape of trajectories out of the aforementioned phase-space region eckmann-ruelle ; kantz-grassberger . This region contains a fractal repeller $`_\chi `$ composed of trajectories which escape neither in future nor in past. The escape rate is characteristic of this fractal repeller and is related to the transport coefficient, leading to the formula (2). The escape-rate formalism has already been applied to the transport property of diffusion gasp-baras  as well as to reaction-diffusion processes claus-gasp ; claus-gasp-vanbei .

The purpose of the present paper is to apply the escape-rate formalism to viscosity. The system we use as a vehicle of our study is a minimal model of viscosity previously analyzed by Bunimovich and Spohn spohn-buni . The minimal models of transport are of special interest because they are the simplest possible models already presenting a positive and finite transport coefficient. It is known that a minimal model should contain only one particle for diffusion, two particles for viscosity, and three particles for heat conduction spohn-buni . For viscosity, we therefore consider here the model composed of two hard disks moving on a torus and undergoing elastic collisions. In a previous paper premier-article , we have described this model and some of its properties for a hexagonal and a square geometries. Our aim is here to compute viscosity thanks to the escape-rate formalism and to show the equivalence with the results of the Green-Kubo and Einstein-Helfand formulas already obtained in the previous paper premier-article .

In the present paper, we compute the viscosity by using the chaotic and fractal properties of the repeller. We use a variant of Eq. (2) in which the difference between positive Lyapunov exponents and KS entropy is given in terms of the Lyapunov exponents and the partial fractal dimensions of the fractal repeller. Indeed, the fractality of the repeller is a corollary of its chaoticity so that its fractal dimensions are related to its KS entropy. Accordingly, the knowledge of the partial dimensions allow us to evaluate the KS entropy.

The paper is organized as follows. In Sec. II, we develop the escape-rate formalism for shear viscosity. In Sec. III, we present the fractal repeller of viscosity in the two-disk model, which we compare with the fractal repeller of diffusion in the Lorentz gas. In this way, we show that a specific fractal repeller is associated with each transport property. Finally, the chaotic and fractal properties are described in Sec. IV where we compute the viscosity coefficient from the positive Lyapunov exponent and the fractal dimension of the repeller. Conclusions are drawn in Sec. V.

## II VISCOSITY IN THE ESCAPE-RATE FORMALISM

### II.1 Helfand moment for viscosity

In the previous paper premier-article , we have shown that the shear viscosity coefficient can be computed with the Einstein-Helfand formula helf

$$\eta =\eta _{xy,xy}=\eta _{yx,yx}=\underset{t\mathrm{}}{lim}\frac{1}{2t}[\stackrel{~}{G}_{yx}(t)\stackrel{~}{G}_{yx}(0)]^2,$$
(3)
where the Helfand moment is defined by

$$\stackrel{~}{G}_{yx}\sqrt{\frac{\beta }{V}}G_{yx}=\sqrt{\frac{\beta }{V}}\left[\underset{a=1}{\overset{N}{}}x_a(t)p_{ay}(t)\underset{a=1}{\overset{N}{}}\underset{s}{}\mathrm{\Delta }x_a^{(s)}p_{ay}^{(s)}\theta (tt_s)\right]$$
(4)
for a system of $`N`$ particles of position $`𝐫_a=(x_a,y_a,\mathrm{})`$ and momentum $`𝐩_a=(p_{ax},p_{ay},\mathrm{})`$ moving in a domain delimited by periodic boundary conditions. As explained in the previous paper premier-article , the particles must satisfy the minimum image convention, which requires the presence of the extra terms in the Helfand moment (4) involving the jumps $`\mathrm{\Delta }x_a^{(s)}`$ of the particles to fulfil the periodic boundary conditions. $`t_s`$ are the times of the jumps. $`p_{ay}^{(s)}=p_{ay}(t_s)`$ is the momentum at the time $`t_s`$ of the jump. In Eq. (3), the average $``$ is taken in the equilibrium microcanonical ensemble for which

$$\beta =\frac{N}{k_\mathrm{B}T(N1)}.$$
(5)
### II.2 First passage problem for viscosity

The central object of the escape-rate formalism is the fractal repeller which is composed of the phase-space trajectories for which the Helfand moment fluctuates forever (in the future and the past) within some interval:

$$\frac{\chi }{2}\stackrel{~}{G}_{yx}+\frac{\chi }{2}.$$
(6)
These trajectories are exceptional because the Helfand moment escape out of this interval for almost all the trajectories. Therefore, the repeller has a vanishing probability measure in the phase space albeit it is typically composed of a nonenumerable set of trajectories. The repeller thus typically forms a fractal in the phase space gasp-dorf ; gasp-book .

We set up a first passage problem of the Helfand moment by introducing absorbing boundaries at $`\stackrel{~}{G}_{yx}=\pm \frac{\chi }{2}`$. These absorbing boundaries in the space of variations of the Helfand moment correspond to equivalent absorbing boundaries in the phase space of the system. In the phase space, the absorbing boundaries delimit a domain which contains the fractal repeller. We consider a statistical ensemble of initial conditions taken inside this domain and we run their trajectories. When a trajectory reaches the absorbing boundaries it escape out of the domain and it thus removed out of the statistical ensemble.

Under the forward time evolution, the remaining trajectories belong to the stable manifolds of the repeller. Under the backward time evolution, the remaining trajectories belong to the unstable manifolds of the repeller. Under both the forward and backward time evolution, the remaining trajectories belong to the repeller itself which is the intersection of its stable and unstable manifolds gasp-dorf . For a typical chaotic dynamics, almost all trajectories escape out of the domain after some time so that the repeller as well as its stable or unstable manifolds are fractal objects.

These fractals can be generated by allowing the escape of trajectories over a long but finite time interval. Over a finite time, there remain a sizable set of trajectories, which progressively reduces to the fractal as the time interval becomes longer and longer. The number of trajectories in the set (or statistical ensemble) decays with time. Typically, the decay is exponential and characterized by the so-called escape rate $`\gamma `$.

The escape rate $`\gamma `$ can be evaluated by solving the first passage problem of the Helfand moment by introducing absorbing boundaries at $`\stackrel{~}{G}_{yx}=\pm \frac{\chi }{2}`$. Indeed, the Einstein-Helfand equation (3) shows that the Helfand moment performs a diffusive-like random walk. Accordingly, the Helfand moment can be considered as a random variable $`g=\stackrel{~}{G}_{yx}`$ for which the probability density $`p(g)`$ obeys to a diffusion-type equation dorf-gasp :

$$\frac{p}{t}=\eta \frac{^2p}{g^2},$$
(7)
where the role of the diffusion coefficient is played by the shear viscosity (3) itself. At the absorbing boundaries, the probability density must satisfy the absorbing boundary conditions:

$$p\left(\frac{\chi }{2}\right)=p\left(+\frac{\chi }{2}\right)=0.$$
(8)
The solution of the diffusion-type equation (7) with the boundary conditions (8) is given by

$$p(g,t)=\underset{n=1}{\overset{\mathrm{}}{}}c_n\mathrm{exp}(\gamma _nt)\mathrm{sin}\left[\frac{\pi n}{\chi }\left(g+\frac{\chi }{2}\right)\right],$$
(9)
with

$$\gamma _n=\eta \left(\frac{\pi n}{\chi }\right)^2,$$
(10)
and where the coefficient $`c_n`$ depend on the initial probability density. The number $`𝒩(t)`$ of trajectories remaining between the absorbing boundaries at the current time $`t`$ is related to the probability density by:

$$𝒩(t)=𝒩_0_{\frac{\chi }{2}}^{+\frac{\chi }{2}}p(g,t)𝑑g𝒩_0\mathrm{exp}(\gamma _1t).$$
(11)
After a long time, the escape is dominated by the smallest decay rate, $`\gamma _1`$, which can therefore be identified with the escape rate $`\gamma `$. In this way, we obtain the escape rate as a function of $`\chi `$:

$$\gamma =\gamma _1=\eta \left(\frac{\pi }{\chi }\right)^2.$$
(12)
This result is obtained by using the diffusion-type equation (7) which is expected to hold over spatial distances larger than the mean free path of the particles. Therefore, the parameter $`\chi `$ must be sufficiently large so that the Helfand moment is in a diffusion regime and Eq. (7) holds.

The shear viscosity coefficient can thus be obtained from the escape rate which depends on the parameter $`\chi `$ of separation between the absorbing boundaries as

$$\eta =\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2\gamma (\chi ).$$
(13)
In the following we call Eq. (13) the escape-transport formula.

### II.3 The chaos-transport formula

At the microscopic level of description, the escape rate is controlled by the fractal repeller $`_\chi `$ which is composed of all the trajectories which satisfy the condition (6) under forward and backward time evolution. The repeller is the support of a natural invariant probability measure. This invariant measure is natural because it is generated by the dynamics and can be approximated by a statistics based on the trajectories remaining within the absorbing boundaries after a long but finite time. The dynamics on the fractal repeller is characterized by positive Lyapunov exponents and a KS entropy, both evaluated with respect to the natural invariant measure of the repeller. If the dynamics is unstable some Lyapunov exponents are positive. If the dynamics is chaotic the KS entropy is positive. On a repeller, the sum of positive Lyapunov exponents differs from the KS entropy and the difference gives the escape rate $`\gamma (\chi )`$ of the repeller $`_\chi `$: eckmann-ruelle ; kantz-grassberger

$$\gamma (\chi )=\left(\underset{\lambda _i>0}{}\lambda _ih_{\mathrm{KS}}\right)__\chi .$$
(14)
If we combine this result from dynamical systems theory with Eq. (13), we obtain the chaos-transport formula (2) for viscosity as originally derived by Dorfman and Gaspard dorf-gasp .

An equivalent formula can be obtained which involves the partial fractal dimensions of the repeller instead of the KS entropy. Indeed, the fractal character of the repeller is a direct consequence of the escape of trajectories so that the KS entropy is no longer equal to the sum of Lyapunov exponents but to

$$h_{\mathrm{KS}}=\underset{\lambda _i>0}{}d_i\lambda _i,$$
(15)
where the coefficients are the partial information dimensions of the repeller associated with each unstable direction of corresponding Lyapunov exponent $`\lambda _i`$ eckmann-ruelle . These partial dimensions satisfy

$$0d_i1,$$
(16)
so that the KS entropy is in general smaller than the sum of positive Lyapunov exponents. Accordingly, the escape rate can be expressed as

$$\gamma (\chi )=\left(\underset{\lambda _i>0}{}c_i\lambda _i\right)__\chi $$
(17)
in terms of the partial codimensions defined as

$$c_i1d_i.$$
(18)
Combining with Eq. (13), the shear viscosity is given by

$$\eta =\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2\left(\underset{\lambda _i>0}{}c_i\lambda _i\right)__\chi .$$
(19)
In the limit $`\chi \mathrm{}`$, the Lyapunov exponents reach their equilibrium values $`\lambda _{i,\mathrm{eq}}`$, while the codimensions vanish typically as $`c_i\chi ^2`$ if transport is normal. If we introduce the quantities

$$a_i\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2c_i|__\chi ,$$
(20)
Eq. (19) provides a decomposition of the viscosity coefficient on the spectrum of Lyapunov exponents such as

$$\eta =\underset{\lambda _{i,\mathrm{eq}}>0}{}a_i\lambda _{i,\mathrm{eq}}.$$
(21)
Typically, the escape is most important in the most unstable direction corresponding to the maximum Lyapunov exponent $`\lambda _1`$. Therefore, the repeller is more depleted in the most unstable direction and the corresponding partial dimension $`d_1`$ is lower than the further ones. This reasoning suggests that a typical behavior is

$$\left(\underset{\lambda _i>0}{}c_i\lambda _i\right)__\chi \left(c_1\lambda _1\right)__\chi ,$$
(22)
for $`\chi \mathrm{}`$ if the maximum Lyapunov exponent $`\lambda _1`$ is well defined.

This is precisely the case in two-degree-of-freedom systems such as the two-disk model where the chaos-transport formula reduces to

$$\eta =\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2\left(c_\mathrm{I}\lambda \right)__\chi ,$$
(23)
where $`\lambda `$ is the unique positive mean Lyapunov exponent and $`c_\mathrm{I}`$ the corresponding codimension which should be understood as the partial information codimension of the unstable manifolds of the fractal repeller given in terms of the partial information dimension by Young’s formula young

$$c_\mathrm{I}=1d_\mathrm{I}=1\frac{h_{\mathrm{KS}}}{\lambda }.$$
(24)
It is known that the partial information dimension of the repeller is well approximated by the partial Hausdorff dimension if the escape rate is small enough and if Ruelle’s topological pressure does not present a discontinuity. This last condition is fulfiled if the system does not undergo a dynamical phase transition, which is the case in the finite-horizon regimes of Sinai’s billiard which controls the dynamics of both the Lorentz gas and the two-disk model gasp-baras . Under these conditions, we can replace the partial information codimension by the partial Hausdorff codimension in the chaos-transport formula and obtain the viscosity as

$$\eta =\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2\left(c_\mathrm{H}\lambda \right)__\chi .$$
(25)
In the limit $`\chi \mathrm{}`$, the Lyapunov exponent converges to its equilibrium value so that we can also write:

$$\eta =\lambda _{\mathrm{eq}}\underset{\chi \mathrm{}}{lim}\left(\frac{\chi }{\pi }\right)^2c_\mathrm{H}(\chi ),$$
(26)
under similar conditions as Eq. (25).

## III Fractal repeller

In this section, our purpose is to display the fractal repeller associated with viscosity in the two-disk model and to compare it with the fractal repeller of diffusion in the Lorentz gas in order to show that they are different and therefore specific to each transport property.

### III.1 Shear viscosity in the two-disk model

The two-disk model has been studied by Bunimovich and Spohn who proved thanks to a central limit theorem that viscosity is positive and finite in this minimal model spohn-buni . In the previous paper premier-article , we have considered the two-disk model in the hexagonal geometry which we shall use in the following. We showed in Ref. premier-article  that the dynamics reduces to a Sinai billiard in the center-of-mass frame and that the Helfand moment (4) with $`N=2`$ is then given by

$$\stackrel{~}{G}_{yx}(t)=\sqrt{\frac{\beta }{V}}\left[x(t)p_y(t)\underset{s}{}\mathrm{\Delta }x^{(s)}p_y^{(s)}\theta (tt_s)\right],$$
(27)
where $`(x,y)`$ are the coordinates of the relative position of both disks and $`(p_x,p_y)`$ the canonically conjugated relative momentum. The jumps happen when the trajectory of Sinai’s billiard crosses the hexagonal boundary. If the trajectory crosses the side of label $`\omega `$ the trajectory is reinjected at the opposite side so that the jump in position is given by the lattice vector $`\mathrm{\Delta }𝐫^{(s)}=𝐜_\omega ^{(s)}`$ corresponding to the side $`\omega `$.

A fractal repeller is defined by considering all the trajectories such that their Helfand moment satisfies the conditions

$$\frac{\chi }{2}\stackrel{~}{G}_{yx}+\frac{\chi }{2},$$
(28)
where the parameter $`\chi `$ should be large enough. The stable manifolds of the fractal repeller can be visualized by plotting the initial conditions of trajectories satisfying the conditions (28) over a long time interval extending forward in time. These initial conditions are taken on the disk of Sinai’s billiard. The initial conditions are specified by the angle $`\theta `$ of the initial position and the angle $`\varphi `$ that the initial velocity makes with a vector which is normal to the disk at the initial position (see Fig. 1). The initial conditions are plotted in the Birkhoff coordinates ($`\theta `$,$`\mathrm{sin}\varphi `$).

Figure 2 depicts such the fractal composed of the stable manifolds of the repeller for viscosity in the two-disk model. We provide evidence that the set is fractal by zooming successively on it in Figs. 3 and 4, where the self-similarity of the repeller clearly appears.

Let us take a section across the repeller in Fig. 2 at $`\theta _0=\pi /4`$. Taking the escape time of the corresponding trajectory, we have obtained the escape-time function depicted in Fig. 5. The time for the trajectory to escape out of the phase-space region corresponding to the interval (28) is infinite if the trajectory belongs to the stable manifold of a trajectory of the repeller. Indeed, this trajectory is then asymptotic to a trajectory which does not escape. Accordingly, the escape-time function has vertical asymptotes on the stable manifolds of the repeller. Since the repeller is fractal the vertical asymptotes are not enumerable, which explains the behavior in Fig. 5.

### III.2 Diffusion in the Lorentz gas

Diffusion of a tracer particle in the hard-disk periodic Lorentz gas has been studied with the escape-rate formalism in Ref. gasp-baras . In this Lorentz gas, the tracer particle undergoes elastic collisions on hard disks forming a triangular lattice. In a unit cell of the lattice, the dynamics also reduces to Sinai’s billiard. The energy of the tracer particle is conserved as well as the phase-space volumes. Sinai and Bunimovich have proved that the dynamics is ergodic and mixing and that the diffusion coefficient is positive and finite in the finite-horizon regime sinai-buni . For diffusion, the associated Helfand moment is simply given by one of the coordinates $`(x,y)`$ of position of the tracer particle dorf-gasp . An escape process is associated with diffusion by setting up a problem of first passage of the tracer particle at some absorbing boundaries. If we consider the $`x`$-coordinates, the tracer particle does not escape as long as the following condition is satisfied:

$$\frac{R}{2}x+\frac{R}{2}.$$
(29)
The absorbing boundary conditions are therefore defined at $`x=\pm \frac{R}{2}`$. With these absorbing boundaries, the system is called an open Lorentz gas gasp-baras .

The trajectories trapped within the interval (29) form a fractal repeller as shown in Ref. gasp-baras . In order to compare with the fractal repeller of viscosity, we can plot the fractal repeller of diffusion in a similar way as here above for viscosity.

Here again, we plot all the initial conditions of trajectories remaining within the interval (29) over a long forward time interval. These initial conditions are plotted in the same Birkhoff coordinates ($`\theta `$,$`\mathrm{sin}\varphi `$) of a disk around the coordinate $`x0`$ in the Lorentz gas. The set of the selected initial conditions approximate the stable manifolds of the fractal repeller. We zoom successively on this fractal in Figs. 7 and 8, which provides evidence of its self-similarity. As a consequence, the repeller is also fractal. The fractal dimension of the repeller is related to the diffusion coefficient of the Lorentz gas and its Lyapunov exponent, as shown in Refs. gasp-nicolis ; gasp-baras .

### III.3 Comparison between diffusion and viscosity

We observe that the two fractal repellers associated respectively with diffusion (see Fig. 6) and viscosity (see Fig. 2) are different. Indeed, although the global structure is similar, the trajectories belonging to the different repellers are not the same.

To convince us of this difference, we take some examples of trajectories. In Fig. 9, we have a periodic trajectory bouncing between two disks in the billiard. This trajectory belongs to the repeller associated with diffusion since the position $`x`$ is bounded and satisfied (29). However, the viscosity Helfand moment of this trajectory does not satisfy the condition (28) so that it does not belong to the repeller of shear viscosity. With Eq. (27), we see that, in one direction, both $`\mathrm{\Delta }x^{(s)}`$ and $`p_y^{(s)}`$ are positive. Therefore, the contribution at this passage is positive for the Helfand moment. In the other direction, both $`\mathrm{\Delta }x^{(s)}`$ and $`p_y^{(s)}`$ are negative but the product $`\mathrm{\Delta }x^{(s)}p_y^{(s)}`$ is also positive. Accordingly, the Helfand moment increases forever on this trajectory which, therefore, does not belong to the repeller associated with shear viscosity.

On the other hand, we can observe the opposite case. Figure 10 depicts an example of trajectory escaping from the interval (29) although its Helfand moment of viscosity remains between the absorbing boundary conditions (28).

The repellers associated with different transport properties are therefore different.

### III.4 Escape rate and viscosity

In this subsection, we show that the shear viscosity can be obtained from the escape rate of the repeller by using the escape-transport formula (13). We consider a sequence of repellers with larger and larger values of the parameter $`\chi `$. The escape rate $`\gamma (\chi )`$ is numerically evaluated for each repeller by computing the decay of the number $`𝒩(t)`$ of trajectories still within the interval (28) at current time and by extracting the escape rate from the exponential decay. The escape rate is observed to behave as $`\gamma (\chi )\chi ^2`$ and the shear viscosity coefficient is then obtained with Eq. (13).

Figures 11 and 12 depict the viscosity directly computed from the escape rate and compared with the values obtained by the Einstein-Helfand formula in Ref. premier-article , respectively in the hexagonal and square geometries. As in the previous paper premier-article , we consider reduced viscosities defined by

$$\eta _{ij,kl}^{}\frac{\eta _{ij,kl}}{2\sqrt{mk_\mathrm{B}T}}.$$
(30)
We observe in Figs. 11 and 12 the excellent agreement between both methods.

## IV Viscosity from the chaotic and fractal properties of the repeller

In this section, we compute the shear viscosity coefficient in terms of the chaotic and fractal properties of the repeller by using the chaos-transport formula (25) which related the viscosity to the Lyapunov exponent and the Hausdorff codimension of the repeller of viscosity.

### IV.1 Lyapunov exponent

In Sinai’s billiard which controls the reduced dynamics of the two-disk model, the elastic collisions between the disks are defocusing. This induces a dynamical instability of the trajectories which is characterized by the Lyapunov exponents. These exponents are the rates of exponential separations between a reference orbit and infinitesimally close orbits. Since the dynamics of Sinai’s billiard is symplectic and volume-preserving in the four-dimensional phase space, the Lyapunov exponents spectrum is ($`+\lambda ,0,0,\lambda `$) so that their sum is vanishing. One of the Lyapunov exponents vanishes because of the absence of exponential separation in the direction of the flow. Another one corresponding to the direction perpendicular to the energy shell equals zero because of energy conservation.

There exists a method to calculate the positive Lyapunov exponent by considering the motion of a front of particles accompanying the reference particle and issued from the same initial position but with different initial velocities sinai . Because the dynamics is defocusing, this front is expanding. Locally on the reference orbit $`\mathrm{\Gamma }_t`$ the front (called the unstable horocycle) is characterized by a curvature $`\kappa _u(\mathrm{\Gamma }_t)`$ or, equivalently, by its radius of curvature $`1/\kappa _u(\mathrm{\Gamma }_t)`$. Thanks to this method explained in detail in Ref. gasp-book ; gasp-baras , we have computed the positive Lyapunov exponent as a function of the density of the system (in the hexagonal and square geometries). The equilibrium values of the Lyapunov exponent are obtained by running a trajectory in Sinai’s billiard without absorbing boundaries and by averaging over a long time interval. The resulting numerical values are depicted in Figs. 13 and 14.

In the chaos-transport formula (25), the Lyapunov exponent has to be evaluated for the trajectories belonging to the fractal repeller. The statistical average is here carried out for the natural invariant probability measure concentrated on the fractal repeller. This invariant measure defines a nonequilibrium state for the motion. As aforementioned, the natural invariant measure is generated by the dynamics itself. Accordingly, the Lyapunov exponent is numerically computed by averaging over a statistical ensemble of trajectories which has not yet escaped after a long but finite time. This ensemble can be as large as wished by increasing the number of initial conditions. In this way, we can calculate the nonequilibrium values of the Lyapunov exponent.

In Table I, we present a comparison between the equilibrium Lyapunov exponent $`\lambda _{\mathrm{eq}}`$ without absorbing boundary conditions (as depicted in Fig. 13 and 14 and the nonequilibrium Lyapunov exponent $`\lambda _{\mathrm{neq}}(\chi )`$ evaluated over a nonequilibrium measure which has the fractal repeller as support. The difference between these exponents is small and of the order of the escape rate, in agreement with the results of Ref. vanbei-latz-dorf  for the disordered Lorentz gas.

Table I. Values of the characteristic quantities of chaos for different densities $`n`$ in the hexagonal system: $`\lambda _{\mathrm{eq}}`$ is the equilibrium Lyapunov exponent for the closed system. The following quantities characterize the fractal repeller for viscosity with $`\chi =60\sqrt{n}`$: $`\lambda _{\mathrm{neq}}`$ is the nonequilibrium Lyapunov exponent of the repeller, $`h_{\mathrm{KS}}`$ its KS entropy, $`\gamma `$ its escape rate, $`c_\mathrm{I}`$ its partial information codimension, and $`c_\mathrm{H}`$ its partial Hausdorff codimension.

| $`n`$ | $`\lambda _{\mathrm{eq}}`$ | $`\lambda _{\mathrm{neq}}`$ | $`\gamma `$ | $`h_{\mathrm{KS}}=\lambda _{\mathrm{neq}}\gamma `$ | $`c_\mathrm{I}=\gamma /\lambda _{\mathrm{neq}}`$ | $`c_\mathrm{H}`$ |
| --- | --- | --- | --- | --- | --- | --- |
| 0.40 | 1.5156 | 1.5163 | 0.0017 | 1.5146 | 0.0011 | 0.0011 |
| 0.50 | 2.3519 | 2.3539 | 0.0023 | 2.3516 | 0.00098 | 0.00092 |
| 0.60 | 3.7258 | 3.7249 | 0.0015 | 3.7234 | 0.00040 | - |

### IV.2 Hausdorff dimension and viscosity

In order to determine the viscosity by the chaos-transport formula (25), we need to determine the partial Hausdorff codimension $`c_\mathrm{H}`$ of the fractal repeller. The corresponding dimension $`d_\mathrm{H}=1c_\mathrm{H}`$ is the Hausdorff dimension of the vertical asymptotes of the escape-time function depicted in Fig. 5. Its values range in the interval $`0d_\mathrm{H}1`$.

The Hausdorff codimension can be obtained using the following numerical algorithm developed by the group of Maryland maryland-group . We consider an ensemble of pairs of trajectories starting from initial conditions $`\varphi _0`$ differing in a value $`ϵ`$. The time taken by the trajectories to escape out of the system is given by the escape-time function in Fig. 5. The pair is said to be uncertain if the trajectories and their Helfand moment present at least one of the following conditions: (i) both trajectories follow paths that differ by the successive passages through the cell boundaries, that is, if we associate to each trajectory a symbolic sequence ($`\omega _1,\omega _2,\mathrm{}`$) which gives the labels of the cell boundaries across which the successive passages occur, and both sequences are different; (ii) one of both trajectories has its Helfand moment which reaches the fixed absorbing boundaries (28) when the Helfand moment of the other one still remains within these limits. If the pair does not present one of these conditions it is called certain. The fraction $`f(ϵ)`$ of uncertain pairs in the initial ensemble is known to behave as the power

$$f(ϵ)ϵ^{c_\mathrm{H}},$$
(31)
giving the Hausdorff codimension as its exponent. Derivations of this result can be found elsewhere ott ; claus-gasp-vanbei ; maryland-group . This method has been already used in Refs. gasp-baras ; claus-gasp ; claus-gasp-vanbei .

We have here applied the Maryland algorithm to obtain the Hausdorff codimension of the fractal repeller of viscosity. Table I compares the partial Hausdorff codimension with the partial information dimension in particular cases. We observe that both codimensions take very close values as expected.

By varying the parameter $`\chi `$, we have obtained the shear viscosity thanks to the chaos-transport formula (25). These values are plotted in Figs. 15 and 16 for the hexagonal and square geometries, respectively. We consider the shear viscosity coefficient $`\eta ^{}`$ in the hexagonal geometry and the element $`\eta _{xy,xy^{}}`$ of the viscosity tensor in the square geometry. The values obtained with the chaos-transport formula (25) are compared with the values obtained by the escape-transport formula (13) and those by the Einstein-Helfand formula (4) obtained in the previous paper premier-article . The agreement between the three formulas is excellent, which confirms the theoretical results.

## V Conclusions

In the present paper, we have applied the escape-rate formalism to the computation of shear viscosity in the two-disk model by Bunimovich and Spohn spohn-buni .

The escape-rate formalism implies the appearance of a fractal repeller associated with viscosity. We have numerically generated the fractal repeller associated with viscosity in this model and we have provided evidence for its fractal character. Using the chaos-transport formula of the escape-rate formalism, we have been able to evaluate the shear viscosity from the positive maximum Lyapunov exponent and the Hausdorff codimension of the fractal repeller of viscosity. The values obtained by using the chaos-transport formula for shear viscosity have been compared with the values obtained by other methods based on the Einstein-Helfand formula, which is equivalent to the Green-Kubo formula as shown in our previous paper premier-article . An excellent agreement has been observed between the different methods. This agreement brings an important support to the escape-rate formalism as a method to establish a connection between the transport properties – here of viscosity – and the underlying microscopic chaotic dynamics. The agreement therefore confirms the theoretical results of the escape-rate formalism gasp-nicolis ; dorf-gasp .

This confirmation opens perspectives for our understanding of the hydrodynamic modes of shear viscosity. Indeed, in the case of diffusion, the fractal character of the repeller has been shown to imply that the nonequilibrium states are described by singular measures defined in the phase space of the system. These results suggest a similar scenario in the case of viscosity that the nonequilibrium states of a sheared fluid should also be described by singular measures at the phase-space level of description.

Our results can be extended to investigate the connection between viscosity and underlying chaos in many-particle systems with a whole spectrum of positive Lyapunov exponents. In many-particle systems, the fractal repeller is characterized by a spectrum of partial fractal dimensions which enter the chaos-transport formula (23). In this way, we can decompose a transport property such as viscosity onto the spectrum of Lyapunov exponents.

The escape-rate formalism also opens perspectives to study viscosity and other transport properties in nanoscopic systems such as nanopipes or atomic and molecular clusters. Indeed, the escape-rate formalism provides a way to define the transport properties already in nanoscopic systems containing a very small number of particles. The escape-rate formalism is particularly appropriate for nanoscopic systems because the transport coefficients can be defined with the escape rate of the Helfand moment out of a finite range of variation. Since the Helfand moment for viscosity can be seen as the center of momenta of the particles its variations are limited to a finite range in a system containing a finite number of particles. In this way, the transport properties could be studied for nanoscopic systems of increasing size in a similar way as the equilibrium thermodynamic properties have been studied in nanoscopic systems as a function of their size.

We hope to report on these issues in future publications.

Acknowledgments. The authors thank Professors J. R. Dorfman and G. Nicolis for support and encouragement in this research, and Dr. I. Claus for helpful discussions. SV thanks the FRIA for financial support. PG thanks the National Fund for Scientific Research (FNRS Belgium) for financial support.


format: text

text:
# Diagonal Antiferromagnetic Easy Axis in Lightly Hole Doped Y1-xCaxBa2Cu3O6

## Abstract

Hole induced changes in the antiferromagnetic structure of a lightly Ca doped Gd:Y1-xCaxBa2Cu3O6 copper oxide single crystal with $`x0.008`$ is investigated by Gd3+ electron spin resonance. Holes do not localize to Ca2+ ions above $`2.5\text{K}`$ since the charge distribution and spin susceptibility next to the Ca2+ are independent of temperature. Both hole doped and pristine crystals are magnetically twinned with an external magnetic field dependent antiferromagnetic domain structure. Unlike the undoped crystal, where the easy magnetic axis is along $`[100]`$ at all temperatures, the easy direction in the hole doped crystal is along the $`[110]`$ diagonal at low temperatures and changes gradually to the $`[100]`$ direction between $`10\text{K}`$ and $`100\text{K}`$. The transition is tentatively attributed to a magnetic anisotropy introduced by hole ordering.

antiferromagnetic domains, antiferromagnetic easy axis, hole localization, stripes

The structure of the magnetic elementary cell is well established in YBa2Cu3O6, the antiferromagnetic (AF) parent compound of a high temperature superconductor, but little is known about how introducing holes affects the magnetic structure. Neutron diffraction found a chessboard like AF order in the basal $`(a,b)`$ plane Tranquada88  which, in high purity crystals, alternates along the $`c`$ direction with the periodicity of the lattice as shown in Fig. 1(a). In YBa2Cu3O6+δ no hole induced change of this structure has been reported up to date. However, in La2-xSrxCuO4, the other well studied AF parent compound for high $`T_\text{c}`$ superconductivity, hole doping changes the magnetic structure in a peculiar way Wakimoto1999 ; Matsuda2000 . At low temperatures a static magnetic modulation appears with a wavelength roughly proportional to the hole concentration and is characterized by a single wavevector running along one of the diagonals of the slightly distorted square CuO2 lattice. It is not yet clear whether this spin density wave is or is not a manifestation of a theoretically predicted Zaanen ; Emery  phase separation into hole-rich and undoped regions. Tranquada et al. Tranquada95  suggested that the stripe like charge and magnetic order in La1.6-xNd0.4SrxCuO4 and related compounds is evidence for such a phase separation. Ando et al. Ando2002  interpreted the doping dependence of the resistivity in terms of segregated conducting stripes of holes. However, to present, electron structural studies Yoshida  do not show the expected quasi 1D Fermi surface planes at low concentrations.

In this work we observe a hole induced change in the AF structure in a single crystal of Gd:Y1-xCaxBa2Cu3O6 with a Ca concentration of $`x0.008`$. Ca2+ substitutes for Y3+ and introduces 0.004 holes per Cu(2) into the AF CuO2 layers of our sample. We show that in the lightly hole doped AF crystal, holes are not localized in the neighborhood of Ca2+ ions down to $`2.5\text{K}`$. We observe AF domains that are static on the time scale of $`10^8\text{s}`$ at temperatures below $`200\text{K}`$. As opposed to La2-xSrxCuO4, YBa2Cu3O6 has a tetragonal crystal structure and is therefore sensitive to any small extra anisotropy. We find that the easy axis of sublattice magnetization rotates by $`45^{}`$ as the temperature is lowered. Unlike the undoped reference that has the structure of Fig. 1(a) at all temperatures, in the doped sample the ground state spin orientation is along the diagonal of the CuO2 square lattice \[Fig. 1(b)\]. The change of the magnetic anisotropy of the bulk by a very small concentration of holes is an indication of an ordered hole structure.

We use 1% of Gd substituting for Y as a weakly perturbing local electron spin resonance (ESR) probe of the magnetic spin susceptibility and charge redistribution in the AF CuO2 sandwich. The Gd3+ ESR $`g`$ factor shifts and zero field splitting (ZFS) parameters are closely related to these quantities. Ref. jan99  details the methods used to precisely determine the $`g`$ factor and the ZFS parameters from the Gd3+ ESR fine structure lines. This technique has proven successful in the study of AF YBa2Cu3O6+δ jan99 ; jan90  and La2CuO4 Rettori  compounds.

Single crystals of Gd:Y1-xCaxBa2Cu3O6 and the reference crystal Gd:YBa2Cu3O6 were grown in BaZrO3 crucibles as described elsewere Erb . The nominal Gd concentration was 1% in agreement with estimates from the intensity of ESR lines from pairs of neighboring Gd3+ ions. The nominal Ca concentration was $`x_{\text{nom}}=0.03`$, but due to the small distribution coefficient of Ca during crystal growth the actual concentration $`x`$ determined by atomic absorption spectroscopy (AAS) is much smaller, between 0.005 and 0.01. A Ca concentration of $`x=0.005`$ is estimated from the Ca satellite ESR intensity. We shall refer to the Ca concentration as $`x=0.008`$.

The ESR spectra testify the high purity of the sample. The powder ESR spectrum of a small amount of a polycrystalline paramagnetic impurity phase (possibly the Y2BaCuO5 “green phase”) and the Gd3+ fine structure (see Fig. 2) of a low concentration (about 0.5% of the Gd3+ ions) of an unidentified Gd site within the single crystal was observed. The X-band ESR spectra of two Ca doped crystals measured at several orientations at $`77\text{K}`$ were the same and one crystal of $`2\times 3\times 0.2\text{mm}`$ was investigated in detail. The ESR spectra of the reference samples reproduced those reported in Ref. jan99  except for a somewhat broader linewidth and a stronger field dependence of the domain structure.

For magnetic field $`B`$ along the $`c`$ axis, a single series of strong fine structure lines (main lines) is observed, corresponding to the majority of Gd3+ ions with only Y3+ ions occupying the rare earth sites in the neighborhood (Fig. 2). These lines, labeled “1” to “7”, are well fitted by simulated spectra using the same ZFS parameters as for the undoped crystal. Several lines with intensities of the order of 1% of the main lines appear at larger amplification. Most of these were observed in undoped samples also and they are due to closely spaced Gd3+-Gd3+ pairs Simon99  (“A1”, “A” and “B1” in Fig. 2). The unidentified Gd site (series “In”) has also been observed in previous studies jan90 .

We assign the satellite lines, “Ca1”, “Ca2”, “Ca6” and “Ca7”, to Gd3+ ions which have a first neighbor rare earth site occupied by a Ca2+ ion. (The inner three lines of the series could not be identified.) This series has not been observed in the undoped crystals. The relative intensities, the temperature dependence of the intensities and the Larmor frequency dependence clearly show that these lines belong to a Gd3+ fine structure series corresponding to a rare earth site with ZFS parameters not very different from the main line series. The strong first neighbor Ca satellite fine structure series of a magnetically aligned superconducting Gd:Y0.85Ca0.15Ba2Cu3O6+δ powder has been easily identified previously Roki-ampere . The shifts of the satellites from the main series in the $`x=0.15`$ sample are similar to those in the $`x=0.008`$ sample studied here.

We first address the question: do the holes localize at low temperatures to the Ca sites or elsewhere? We compare the temperature dependence of the charge distribution and spin susceptibility at sites next to and distant from Ca2+ ions. Gd3+ ZFS parameters are highly sensitive to the charge distribution. Doping Gd:YBa2Cu3O6+δ with oxygen from $`\delta =0`$ to $`\delta =1`$ changes the largest ZFS parameter, $`b_2^0`$, by 40% Roki92  and changes of this order are expected for the first neighbor Gd3+ if holes were to localize to the Ca2+ ions within a few lattice constants. Hole localization would then shift the “Can” satellite lines on the order of their distances from the central (“4”) line. This is not the case. The observed small, few %, temperature dependence of the ZFS parameters (Fig. 3, squares) is due to thermal expansion of the lattice since the temperature variation is the same at the first neighbors of Ca2+ and at sites distant from the Ca2+ ions. The extra ZFS caused by the Ca2+ neighbor (Fig. 3, triangles), measured by the difference between the “Ca2” (“Ca6”) and “2” (“6”) line positions, is temperature independent within $`\pm 2`$% from 2.5 to $`140\text{K}`$. Thus holes do not localize at Ca2+ ions (at least above $`2.5\text{K}`$) since the charge redistribution at Cu(2) and O(2) lattice atoms surrounding the Ca2+ ions would inevitably change the ZFS. Also, the $`c`$ axis susceptibilities (measured by the $`g`$ factors) near the Ca site and in the bulk are the same within $`\pm 10`$% in the temperature range of 2.5 to $`140\text{K}`$ (Fig. 3) and this reinforces the conclusion that localization of holes does not take place at Ca sites. The equal susceptibilities at the satellite and main sites exclude the possibility that charges are already localized to Ca2+ ions at high temperatures.

As shown below, even a light hole doping has a dramatic effect on the magnetic structure. The temperature independent AF domain structure Burlet ; jan99  of pristine YBa2Cu3O6 was confirmed for the reference sample. The easy direction of the sublattice magnetization $`𝑴_\text{s}`$ is along $`[100]_t`$ (any of the equivalent $`[100]`$ or $`[010]`$ directions in the tetragonal crystal) in the undoped antiferromagnet. In zero applied field, single crystals are magnetically twinned with equal amounts of domains in the two orthogonal directions, $`[100]_t`$, of the tetragonal lattice. The domain structure is field dependent, and applying a magnetic field of a few T along a $`[100]_t`$ direction turns all domains perpendicular to the applied field. In the present reference sample about 80% of domains had $`𝑴_\text{s}`$ perpendicular to the field in a field of $`1\text{T}`$ applied along $`[100]_t`$.

Nothing is known from experiments about the direction, width, concentration or origin of the domain walls separating the orthogonal domains. Domain walls in the $`(a,b)`$ plane that separate domains along $`c`$ were suggested to explain data in undoped crystals jan99 . Conducting stripes running along $`a`$ and separating domains within the $`(a,b)`$ plane were suggested by Niedermayer et al. Niedermayer  and Ando et al. Ando99  to explain NMR and magnetoresistance data.

The structure of the ESR lines reflects the distribution of domain orientations. There are two contributions to the Gd3+ resonance frequency that depend on the orientation of $`𝑴_\text{s}`$: one formally described by a “$`g`$ shift” and the other by an orthorhombic ZFS parameter $`b_2^2`$. Most probably, both are due to the exchange interaction $`J`$ between the CuO2 planes and Gd3+ ions and are proportional to the anisotropic AF spin susceptibility. $`b_2^2`$ may also arise from magnetostriction, but the experiments described here suggest that a term second order in $`J`$, overlooked in Refs. jan99  and Rettori , is dominant.

The AF domain structure is readily apparent in the Gd3+ ESR spectra with magnetic field applied in the $`(a,b)`$ plane. Fig. 4 displays one component of the fine structure of the $`9\text{GHz}`$ ESR spectra at several temperatures for magnetic fields oriented along $`[100]_t`$ (i.e., magnetic field $`𝑩`$ oriented along $`\varphi _𝑩=0^{}`$ or $`90^{}`$) and $`[110]_t`$ ($`\varphi _𝑩=45^{}`$ or $`135^{}`$). At high temperatures the spectra of the Ca doped and undoped YBa2Cu3O6 crystals are similar: the main series are split by the same amount and the line widths are comparable. In contrast to the reference, the ESR spectrum of Gd:Y1-xCaxBa2Cu3O6 is temperature dependent. A gradual magnetic phase transition is observed around $`T_m=40\text{K}`$. Spectra at $`TT_m`$ and $`TT_m`$ are qualitatively different. The transition is broad and non trivial changes with temperature are still observed above $`85\text{K}`$ and below $`15\text{K}`$. The domain structure and the transition depend on magnetic field. An account of the transition between 0.2 and $`5.6\text{T}`$ will be presented elsewhere.

At lower temperatures the magnetic easy axis is rotated by $`\pi /4`$. Below about $`20\text{K}`$, the $`\varphi _𝑩`$ dependence of the lines is shifted by $`\pi /4`$, i.e., the lines in the $`𝑩[100]_t`$ spectrum are similar to the corresponding ones at high temperature in the $`𝑩[110]_t`$ spectrum. We denote the high and low temperature orders as AFI(0) and AFI($`\pi /4`$), respectively. The probable cell for AFI($`\pi /4`$) is shown in Fig. 1(b), but other, larger magnetic cells are also possible. The high temperature spectra are well modeled by the same spin Hamiltonian with the same ZFS parameters as for the undoped case, and this is still true at low temperatures except for a $`\pi /4`$ rotation of the orthorhombic term $`b_2^2`$.

At high temperatures the system is purely AFI(0) within experimental uncertainties while it is almost purely AFI($`\pi /4`$) at the lowest temperatures. During the transition the system is inhomogeneous but it is not possible to tell from an analysis of the ESR spectra whether during the phase transition domains rotate continuously or a mixture of the two phases occurs. The lack of hysteresis with cycling the magnetic field and temperature makes a continuous transition more probable.

The AFI(0) to AFI($`\pi /4`$) phase transition observed here in a high purity, lightly Ca doped crystal has no relation with the AFI(0) to AFII phase transition observed in strongly doped crystals where dopants give rise to magnetic moments in the Cu(1) layer Kadowaki ; Brecht . In the AFII phase every second bilayer is rotated by $`\pi `$ and the lattice is doubled along $`c`$. This could be easily distinguished from the AFI($`\pi /4`$) state by ESR.

We believe the transition is related to hole localization in spite of the low hole concentration $`p=0.004`$. The sample is of high purity and extrinsic effects are unlikely. The broad temperature range of the magnetic transition (i.e., a continuous change in the distribution of the orientations of $`𝑴_\text{s}`$ between 100 and $`6\text{K}`$) is consistent with the idea that the magnetic transition is related to the localization of holes. The gradual localization of holes is well documented in Ca doped YBa2Cu3O6 Niedermayer . At low temperatures the sublattice magnetization of the bulk becomes the same in hole doped systems as in undoped ones, and for Ca concentrations $`x=0.03`$ to 0.06 a broad temperature range of magnetic fluctuations and a magnetic transition to a static order in the range of 10–$`20\text{K}`$ is observed.

Our findings suggest that in Y1-xCaxBa2Cu3O6 in analogy to lightly doped La2-xSrxCuO4, holes localize into an ordered structure which induces a diagonal spin density wave with spin direction coupled to the propagation vector. The amplitude of the spin density wave is small and the change of the easy axis direction is not accompanied by an observable change of the magnitude of the magnetic susceptibility. However, localized holes change both the direction and the magnitude of the magnetic crystalline anisotropy, since at low temperatures the domain structure persists to high magnetic fields. In the Ca doped sample AF domains are well observed at $`5.4\text{T}`$ while they are hardly discernible in the undoped reference at $`2.7\text{T}`$.

It is probable that the easy axis in pristine YBa2Cu3O6 is determined by some extrinsic effects, e.g., residual Cu-O chain fragments. This could explain why no change in the easy direction with temperature was observed in oxygen doped samples Casalta . The difference between hole doping by Ca and oxygen is that oxygen chain fragments introduce a local orthorhombic distortion of the structure which at higher concentrations may fix the direction of magnetic crystalline anisotropy at all temperatures.

In conclusion the main result is a new magnetic phase transition in lightly Ca doped YBa2Cu3O6 in which the easy direction of the AF sublattice magnetization is rotated by $`\pi /4`$ in the $`(a,b)`$ plane. We find that holes do not localize to the Ca dopants suggesting that the transition is related to the condensation of holes into an ordered structure. In view of the similarities of the phase diagrams Niedermayer , the diagonal magnetic order in Y1-xCaxBa2Cu3O6 and the diagonal incommensurate spin modulation observed in doped LaxSr2-xCuO4 may have a common origin.

###### Acknowledgements.

We are indebted to Patrik Fazekas and Peter Littlewood for illuminating discussions and to W. Wendl (Universität Karlsruhe) for the AAS analysis. Support of the Hungarian state grants OTKA TS040878, T029150, T043255 and the European Infrastructure Network, SENTINEL, are acknowledged. The work in Lausanne was partly supported by the NCCR research pool “MaNEP” of the Swiss National Science Foundation.


text:
# MEASUREMENT OF THE ATMOSPHERIC MUON SPECTRUM FROM 20 TO 2000 GeV

## 1 Introduction

Atmospheric muons are amongst the final products of primary cosmic ray induced air shower cascades. A precise measurement of the ground level muon flux can therefore be used to test the understanding of the primary cosmic ray flux and the hadronic interactions involved in the production of the muons’ parent mesons. Moreover, it provides a crucial test for theoretical neutrino flux calculations, because each muon is produced with an accompanying muon neutrino.
Here a new measurement of the atmospheric muon flux is presented using the precise muon spectrometer of the L3 detector  located at the LEP accelerator at CERN, Geneva. During the analysis special attention was given to the precise determination of all relevant detector and environmental parameters needed to convert the raw data distributions to an absolute surface level flux. Due to the large amount of available statistics, extensive studies of the residual systematic uncertainties were possible.

## 2 Experimental Setup

The momentum distribution of atmospheric muons is measured with the extended L3 setup known as L3+C . It is located 450 m above sea level and shielded from the hadronic and electromagnetic air shower components by a 30 m thick molasse overburden. The muon momentum is measured with the L3 muon spectrometer, which is situated inside a large magnetic volume of 1000 m3 at a field of 0.5 T. A 202 m2 scintillator array was installed on top of the detector to record the arrival time of the muons. With this arrangement a relative momentum resolution ranging from 2.2% at 20 GeV to 52% at 2000 GeV is achieved. Being equipped with a trigger and data acquisition system independent of the normal L3 data taking, L3+C recorded 1.2$`10^{10}`$ atmospheric muon triggers during its operation in the years 1999 and 2000.

## 3 Analysis

In total 2$`10^7`$ high quality events are used in the muon spectrum analysis. The geometrical acceptance of the detector as well as the energy loss in the molasse overburden is calculated with a simulation of the L3+C setup and its surroundings. The detector efficiencies are determined from the data itself as a function of time, charge, momentum and zenith angle. The raw event distributions at the detector level are deconvoluted taking into account the detector resolution and the stochastic energy loss in the molasse. The systematic uncertainties of the muon flux and charge ratio measurement are studied by investigating the stability of the results with time and azimuth angle and under the variation of the selection cuts.

## 4 Results

The L3+C vertical muon spectrum is shown in figure 1(a) along with previous measurements  providing an independent absolute normalization. The best precision is achieved around 100 GeV, where the total error amounts to 2.6%. At lower energies uncertainties of the molasse overburden get important, whereas at high energies the statistical error dominates. The measured vertical charge ratio is compared to an average  of previous measurements in figure 1(b). It is worthwhile pointing out, that the total error of one zenith angle bin from this experiment is compatible with the error from the average of all previous experiments. Both, muon flux and charge ratio were measured as function of the zenith angle from 0 to 58 degrees. Their relative angular dependence is found in good agreement with the prediction from an air shower simulation  of atmospheric muons.


text:
# I. Introduction

## I. Introduction

Radial momentum operators naturally arise in quantum mechanics whenever the Hamiltonian contains a radial potential term. In spherical coordinates

$$\mathrm{}^2=r^1_r^2r+L^2/r^2$$

where $`𝐋`$ is the angular momentum operator. The squared radial momentum operator

$$p_r^2r^1_r^2r$$
$`(1)`$
is positive and Hermitian with respect to the usual scalar product space in spherical coordinates

$$\varphi __1|\varphi __2_\mathrm{\Omega }_0^{\mathrm{}}\varphi __1^{}\varphi __2r^2𝑑r𝑑\mathrm{\Omega }.$$
$`(2)`$
The bounded eigenfunctions of $`p_r^2`$ satisfying

$$p_r^2\varphi =k^2\varphi $$
$`(3)`$
are

$$\varphi _k=\mathrm{sin}(kr)/r.$$
$`(4)`$
The singular functions $`[\mathrm{cos}(kr)/r]`$ are excluded as solutions of $`\mathrm{}^2\varphi =k^2\varphi `$ because $`\mathrm{}^2[\mathrm{cos}(kr)/r]`$ results in a $`\delta `$ function at the origin as was pointed out by Dirac (p154 of ), so one imposes the boundary condition

$$r\varphi =0\text{for }r=0.$$
$`(5)`$
The necessity of this boundary condition arises from the singular nature of the spherical coordinate system at the origin.

Various authors have pointed out that the operator

$$\stackrel{~}{p}_rir^1_rr=\frac{1}{2}\left[\frac{𝐫}{r}(i\mathrm{})+(i\mathrm{})\frac{𝐫}{r}\right]$$

is not Hermitian despite its symmetric structure: the underlying reason being that the radial integration within the scalar product (2) is over the half axis instead of the whole axis. Supposing for the moment that $`\stackrel{~}{p}_r`$ is Hermitian, then there exists the unitary operator $`[\mathrm{exp}(ia\stackrel{~}{p}_r)]`$ generating the transformation $`rra,`$ so the spectrum of $`r`$ must extend from $`\mathrm{}`$ to $`\mathrm{}`$ which is clearly false as $`r`$ is bounded below . And the displacement $`rra`$ is clearly unphysical as it would tear a hole at the origin. For an extensive discussion of the non-Hermiticity of $`\stackrel{~}{p}_r,`$ see .

For a positive Hermitian operator such as $`p_r^2`$ there exists a unique square root Hermitian operator which we will call $`p_r^+`$ such that $`p_r^+p_r^+=p_r^2`$ and $`p_r^+`$ is positive. To simplify matters we will eliminate the measure $`r^2`$ in the scalar product (1) by making the substitution

$$\varphi (r)=\chi (r)/r,$$

then the eigenvalue equation (3) is simply

$$_r^2\chi =k^2\chi $$

and the scalar product space is $`L^2(0,\mathrm{}):`$

$$\chi __1|\chi __2_r_0^{\mathrm{}}\chi __1^{}\chi __2𝑑r$$
$`(6)`$
with the boundary condition corresponding to (5)

$$\chi (0)=0.$$
$`(7)`$
Once we have found the square root operator of $`_r^2`$ which we will call $`z_r^+,`$ which is Hermitian with respect to the space (6), then

$$p_r^+=(\frac{1}{r}z_r^+r).$$
$`(8)`$
We will show that the required operator $`z_r^+`$ is

$$\begin{array}{ccc}\hfill z_r^+& =_e_r=_r_o\hfill & (9)\hfill \end{array}$$

where $`_e,_o`$ are the Hilbert transforms of even, odd functions defined by

$$\begin{array}{cc}\hfill _ef(r)& \frac{2r}{\pi }_0^{\mathrm{}}\frac{f(t)}{r^2t^2}𝑑t\hfill \\ \hfill _of(r)& \frac{2}{\pi }_0^{\mathrm{}}\frac{tf(t)}{r^2t^2}𝑑t,\hfill \end{array}$$

the integral being a Cauchy principle value. It is known that $`_e`$ is the Fourier sine transform of the Fourier cosine transform, while $`_o`$ is the Fourier cosine transform of the Fourier sine transform, i.e.

$$\begin{array}{cc}\hfill _e& _s_c\hfill \\ \hfill _o& _c_s\hfill \end{array}$$
$`(10)`$
where

$$\begin{array}{cc}\hfill _sf(r)& \sqrt{\frac{2}{\pi }}_0^{\mathrm{}}f(t)\mathrm{sin}(rt)𝑑t\hfill \\ \hfill _cf(r)& \sqrt{\frac{2}{\pi }}_0^{\mathrm{}}f(t)\mathrm{cos}(rt)𝑑t.\hfill \end{array}$$

Noting the operator identities

$$\begin{array}{cc}\hfill _s_r& =r_c\hfill \\ \hfill _c_x& =+r_s\hfill \end{array}$$
$`(11)`$
then substituting (10) and (11) into (9) we obtain

$$\begin{array}{ccc}\hfill z_r^+& _e_r=_s_c_r=_sr_s.\hfill & (12)\hfill \end{array}$$

And

$$\begin{array}{ccc}\hfill (z_r^+)^2& =(_sr_s)(_sr_s)=_sr^2_s=_r^2\hfill & (13)\hfill \end{array}$$

which demonstrates the required property $`(z_r^+)^2=_r^2.`$ For (12) we have used the identity $`_s_s=\mathrm{\hspace{0.17em}1}=_c_c.`$ The result (12) can also be obtained more directly from (9) using the identity $`_e_o=_o_e=1.`$

Eigenfunctions of $`z_r^+`$  It can be verified that the eigenfunction $`\mathrm{sin}(kr)`$ of $`(_r^2)`$ is also an eigenfunction of $`z_r^+,`$ as

$$\begin{array}{ccc}\hfill z_r^+\mathrm{sin}(kr)& _e_r\mathrm{sin}(kr)=k_e\mathrm{cos}(kr)=k\mathrm{sin}(kr),k>0.\hfill & \end{array}$$

Note that $`\mathrm{cos}(kr)`$ is not an eigenfunction of $`z_r^+,`$ in fact

$$\begin{array}{ccc}\hfill z_r^+\mathrm{cos}(\lambda r)& =k_e\mathrm{sin}(kr)=\frac{2k}{\pi }[\mathrm{sin}(kr)\text{Ci }(kr)\mathrm{cos}(kr)\text{Si }(kr)].\hfill & \end{array}$$

$`z_r^+`$ is a positive operator  The fact that $`z_r^+`$ is a positive Hermitian operator can be seen from the result (12) which is $`z_r^+=_sr_s.`$ The operator $`_s`$ is Hermitian with respect to the $`L^2(0,\mathrm{})`$ space (6), which can be shown by writing out $`\psi |_s\chi .`$ then changing the order of integration, i.e.

$$\begin{array}{ccc}\hfill \psi |_s\chi & _0^{\mathrm{}}\psi ^{}(r)\left[_0^{\mathrm{}}\mathrm{sin}(rt)\chi (t)𝑑t\right]𝑑r\hfill & \\ & =_0^{\mathrm{}}\chi (t)\left[_0^{\mathrm{}}\psi ^{}(r)\mathrm{sin}(rt)𝑑r\right]𝑑t=_s\psi |\chi .\hfill & (14)\hfill \end{array}$$

Then using (14)

$$\begin{array}{ccc}\hfill \chi |z_r^+\chi & =\chi |_sr_s\chi =\sqrt{r}_s\chi |\sqrt{r}_s\chi \mathrm{\hspace{0.17em}0}.\hfill & (15)\hfill \end{array}$$

The inverse of $`z_r^+`$  Every positive operator has an inverse which is also positive. From (12) we see that

$$\begin{array}{ccc}\hfill (z_r^+)^1& =_s\frac{1}{r}_s.\hfill & (16)\hfill \end{array}$$

which is clearly positive by similar reasoning to (15). Writing out (16)

$$\begin{array}{ccc}\hfill (z_r^+)^1\chi (r)& _0^{\mathrm{}}𝑑t\mathrm{sin}(rt)\frac{1}{t}_0^{\mathrm{}}\mathrm{sin}(ut)\chi (u)𝑑u\hfill & \\ \multicolumn{3}{c}{andcarryingouttheintegrationwithrespectto\text{ t }weobtain}\\ \hfill (z_r^+)^1\chi (r)& =\frac{1}{2}_0^{\mathrm{}}\mathrm{log}\left|\frac{r+u}{ru}\right|\chi (u)𝑑u.\hfill & \end{array}$$

$`(z_r^+)^1`$ in terms of the fractional integral operators  Erdelyi  and Kober  established relations between the Hankel transforms (of which $`_s,_c`$ are special cases) and the so-called fractional integral operators (sometimes called Erdelyi-Kober operators after their originators). In particular

$$\begin{array}{ccc}\hfill (z_r^+)^1& =_s\frac{1}{r}_s=r𝒦\hfill & (17)\hfill \end{array}$$

where the Erdelyi-Kober operators $`,𝒦`$ are defined by

$$\begin{array}{cc}& f(r)\frac{2}{\sqrt{\pi }}_0^1u(1u^2)^{1/2}f(ur)𝑑u\hfill \\ & 𝒦f(r)\frac{2}{\sqrt{\pi }}_1^{\mathrm{}}(u^21)^{1/2}f(ur)𝑑u.\hfill \end{array}$$
$`(18)`$
Due to the Hermitian conjugate relations

$$(r𝒦)^{\mathrm{}}=r(r)^{\mathrm{}}=r𝒦$$

then

$$\begin{array}{ccc}\hfill \psi |(z_r^+)^1\chi & =\psi |r𝒦\chi =r\psi |\chi \hfill & (19)\hfill \end{array}$$

which is another way of showing that $`(z_r^+)^1`$ is positive. If we insert the eigenfunctions of $`z_r^+`$ into (19), with $`\chi =\mathrm{sin}(kr),\psi =\mathrm{sin}(k^{}r),`$ and as

$$\mathrm{sin}(kr)=\sqrt{\pi }J_1(kr)$$

we arrive at

$$\begin{array}{ccc}\hfill \mathrm{sin}(k^{}r)|(z_r^+)^1\mathrm{sin}(kr)& =\pi rJ_1(k^{}r)|J_1(kr)=\frac{\pi }{k}\delta (kk^{})\hfill & \end{array}$$

which is an interesting relation between the orthogonality of sine functions and Bessel functions. The operators $`,𝒦`$ are invertible, for a succinct introduction to these operators see chapter 2 of . Rooney  used the operator identities

$$\begin{array}{cc}\hfill _e& =𝒦^1r\frac{1}{r}_o=r^1\frac{1}{r}𝒦\hfill \end{array}$$

to investigate the ranges of $`_e,_o.`$

Conclusion  While there has been discussion in the literature of the $`\stackrel{~}{p}_r`$ operator and its non-Hermiticity, the $`p_r^+`$ operator

$$p_r^+=\frac{1}{r}z_r^+r=\frac{1}{r}_e_rr$$

appears to be little known. Unfortunately the $`p_r^+`$ operator has no simple commutation relation with $`r.`$ A possible application of this operator is in light-cone quantum theory, where the Hamiltonian contains an inverse radial momentum operator - see  which addresses this problem and the references therein, or the original work by Dirac .

## References

 Dirac P. A. M., The Principles of Quantum Mechanics (OUP, Oxford, 1947) 3rd ed.

 Levin O. and Peres A., “Some oddities of lightcone dynamics,” J. Phys. A 27, L143-145 (1994)

 Liboff R. L. et al, “On the radial momentum operator,” Am. J. Phys. 41, 976-980 (1973)

 Erdelyi A., “On fractional integration and its application to the theory of Hankel transforms,” Quart. J. of Math. Oxford 11, 293-303 (1940)

 Kober H., “On fractional integrals and derivatives,” Quart. J. of Math. Oxford 11, 193-211 (1940)

 Sneddon I. N., Mixed boundary problems in potential theory (North Holland, Amsterdam, 1966)

 Rooney P. G., “On the ranges of certain fractional integrals,” Can. J. Math. 24, 1198-1216 (1972)

 Dirac P. A. M., “Forms of relativistic dynamics,” Rev. Mod. Phys. 21, 392-399 (1949)


text:
# MEM Analysis of Glueball Correlators at 𝑇>0

## Abstract

Maximum Entropy Method (MEM) is applied to glueball correlators at finite temperature constructed by using 5,500 – 9,900 gauge field configurations generated in SU(3) quenched lattice QCD with the lattice parameter $`\beta _{\mathrm{lat}}=2N_c/g^2=6.25`$ and the renormalized anisotropy $`a_s/a_t=4`$. The results support the thermal width broadening of the $`0^{++}`$ glueball near the critical temperature.

It is believed that, in the neighborhood of the critical temperature $`T_c`$, the QCD vacuum begins to change its structure such as the reduction of the string tension, the partial restoration of spontaneous chiral-symmetry breaking and so on. Since the hadrons are composites of quarks and gluons, whose interactions depend on the properties of the QCD vacuum, some of the hadrons are expected to change their structure drastically near $`T_c`$. For instance, effective model studies suggest the pole-mass reductions of charmonium and light $`q\overline{q}`$ mesons as consequences of the reduction of the string tension and the partial restoration of spontaneous breaking of chiral symmetry, respectively . In fact, these changes are considered as important precritical phenomena of QCD phase transition in RHIC QGP experiments, and corresponding lattice QCD calculations have been performed at quenched level . In this paper, we study the thermal $`0^{++}`$ glueball, whose mass reduction seems natural in terms of the large difference between the glueball mass $`m_\mathrm{G}`$ 1.5GeV and $`T_c`$ 280MeV in quenched QCD .

To study the pole mass of a hadron in lattice QCD, one has first to construct a temporal correlator as $`G(\tau )=\varphi (\tau )\varphi (0)`$, and then resorts to its spectral representation as

$$G(\tau )=_0^{\mathrm{}}𝑑\omega K(\tau ,\omega )A(\omega ),$$
(1)
where $`K(\tau ,\omega )\frac{\mathrm{cosh}(\omega (\beta /2\tau ))}{\mathrm{sinh}(\beta \omega /2)}`$, $`\beta 1/T`$, and $`A(\omega )`$ is the spectral function with its spatial momentum projected to zero, i.e., $`A(\omega )A(\omega ,\stackrel{}{p}=\stackrel{}{0})`$. Each peak position of $`A(\omega )`$ provides us with a pole mass of a hadron at $`T>0`$, which can be observed in high-energy experiments.

The most popular way to extract the spectral function $`A(\omega )`$ from the constructed correlator is to adopt an ansatz for $`A(\omega )`$, and then to perform the fitting analysis. In the case of glueballs, we first investigated their pole masses using the simple narrow-peak ansatz with $`A(\omega )=C\left\{\delta (\omega m)\delta (\omega +m)\right\}`$, where $`C`$ and $`m`$ are the strength and the temperature-dependent pole mass, respectively .

We then proceeded to an advanced analysis adopting the Breit-Wigner ansatz, a more general and sophisticated ansatz, for the shape of the spectral function as

$`A(\omega )=C\left\{\delta _\mathrm{\Gamma }(\omega \omega _0)\delta _\mathrm{\Gamma }(\omega +\omega _0)\right\}`$ (2)
expressed with the peak function of

$`\delta _\mathrm{\Gamma }(\omega ){\displaystyle \frac{1}{\pi }}{\displaystyle \frac{\mathrm{\Gamma }}{\omega ^2+\mathrm{\Gamma }^2}}.`$ (3)
Here, $`C`$, $`\omega _0`$ and $`\mathrm{\Gamma }`$ are the strength, the peak center corresponding to the pole mass, and the thermal width of the glueball peak, respectively. (The thermal width should not be confused with the decay width. It is important to keep in mind that, at $`T>0`$, the thermal width is generated even for stable particles through the interaction with the thermally-excited particles.)

The Breit-Wigner analysis indicates that the main thermal effect for the $`0^{++}`$ glueball appears as the thermal width broadening  of about 300 MeV with a reduction in the peak center of about 100 MeV. This Breit-Wigner ansatz for the lowest-lying peak in the spectral function seems natural and solid at relatively low temperature. However, near and above the critical temperature, more complicated structure may appear in the spectral function.

In this point, Maximum Entropy Method (MEM) is quite attractive, because it provides us with a numerical procedure to reconstruct the spectral function $`A(\omega )`$ directly from lattice QCD Monte Carlo calculations . In this paper, we present our preliminary results on the MEM-reconstructed $`0^{++}`$ glueball spectrum at $`T>0`$.

We first consider the smearing method. Use of the extended operator is one of the most useful techniques to enhance the low-energy spectra. However, it has a disadvantage that it may make the physical interpretation of resulting peak less trivial. It can create an unphysical bump structure even in the spectral functions of non-interacting particles . For instance, we consider the extreme case for glueballs. The naive continuum limit of the $`0^{++}`$ glueball operator is given as $`\varphi G_{ij}^aG_{ij}^a`$. After the smearing, it is spatially extended according to the Gaussian distribution as

$$\varphi ^{}\begin{array}{c}\frac{d^3yd^3z}{(2\pi )^{3/2}\rho ^3}\mathrm{exp}\left(\frac{(\stackrel{}{y}\stackrel{}{z})^2}{2\rho ^2}\right)G_{ij}^a(\stackrel{}{y})G_{ij}^a(\stackrel{}{z}),\hfill \end{array}$$
(4)
where $`\rho `$ controls the size of the spatial extension . We adopt here the Coulomb gauge to obtain the continuum expression. We calculate the corresponding two-point functions perturbatively up to $`O(\alpha _s^0)`$, and extract their spectral functions. While the non-smeared spectral function is given as $`A(\omega )\omega ^4`$, the smeared spectral function acquires an additional Gaussian factor as

$$A(\omega )\omega ^4\mathrm{exp}\left\{(\omega \rho )^2/4\right\},$$
(5)
which possesses an unphysical bump at $`\omega =2\sqrt{2}/\rho `$. The smearing method suppresses the overlap to the higher spectral components, and, as a consequence, this factor appears.

The actual low-lying glueball is not a perturbative system but a definite bound state/resonance in quenched QCD below $`T_c`$. Note that, in order to study the mass and the width of a difinite resonance, the problem of unphysical bump is not serious, because the pole position in the complex $`\omega `$ plane is unaffected by the smearing . In this paper, we use the smearing method to enhance the lowest-lying $`0^{++}`$ glueball peak, which is negligibly small in the non-smeared spectral function.

We use the SU(3) anisotropic lattice plaquette action  as

$$S_\mathrm{G}=\begin{array}{c}\frac{\beta _{\mathrm{lat}}}{N_c}\frac{1}{\gamma _\mathrm{G}}\underset{s,i<j3}{}\text{Re}\text{Tr}\left\{1P_{ij}(s)\right\}\hfill \\ +\frac{\beta _{\mathrm{lat}}}{N_c}\gamma _\mathrm{G}\underset{s,i,j3}{}\text{Re}\text{Tr}\left\{1P_{i4}(s)\right\},\hfill \end{array}$$
(6)
where $`P_{\mu \nu }(s)\text{SU}(3)`$ denotes the plaquette operator in the $`\mu `$-$`\nu `$-plane. The lattice parameter and the bare anisotropic parameter are fixed as $`\beta _{\mathrm{lat}}2N_c/g^2=6.25`$ and $`\gamma _\mathrm{G}=3.2552`$, respectively, so as to reproduce the renormalized anisotropy as $`a_s/a_t=4`$. The scale unit is introduced from the on-axis data of the static inter-quark potential with the string tension $`\sqrt{\sigma }=440`$ MeV. The resulting lattice spacings are given as $`a_t^1=9.365(66)`$ GeV and $`a_s^1=2.341(16)`$ GeV. The critical temperature is estimated as $`T_c280`$ MeV from the susceptibility of the Polyakov loop. We generate 5,500–9,900 gauge configurations to construct the glueball correlators, where statistical data are divided into bins of the size 100 to reduce possible auto-correlations near the critical temperature. We adopt an appropriate smearing to enhance the low-energy spectrum.

To reconstruct the spectral function $`A(\omega )`$, we apply MEM to the appropriately smeared glueball correlator normalized according to $`G(\tau =0)=1`$. As the practical numerical procedure, we in principle follow Ref. . We adopt the following Shannon-Jaynes entropy as

$$S_0^{\mathrm{}}\left[A(\omega )m(\omega )A(\omega )\mathrm{log}\left(\frac{A(\omega )}{m(\omega )}\right)\right],$$
(7)
where $`m(\omega )`$ is a real and positive function referred to as the default model function. $`m(\omega )`$ is required to mimic the asymptotic behavior of $`A(\omega )`$ as $`\omega \mathrm{}`$. As $`m(\omega )`$, we adopt Eq.(5), i.e., the perturbative spectral function up to $`O(\alpha _s^0)`$ as

$$m(\omega )=N\omega ^4\mathrm{exp}\left\{(\omega \rho )^2/4\right\},$$
(8)
where the normalization factor $`N`$ is determined so as to mimic $`G(\tau =0)=1`$, i.e.,

$$1=_0^{\mathrm{}}𝑑\omega K(\tau =0,\omega )m(\omega ).$$
(9)
In Fig. 1, we show the reconstructed spectral functions of the lowest $`0^{++}`$ glueball from the temporal correlations at $`T=130,253,275`$ MeV. Since the error bar estimated by following Ref.  appears to be unreasonably small, we do not put it in Fig. 1 to avoid unnecessary confusion. For a reasonable estimate, it seems necessary to use the jackknife error estimate . In Fig. 1, we see the tendency that the peak becomes broader with increasing temperature.

To summarize, we have applied Maximum Entropy Method (MEM) to appropriately smeared glueball correlators constructed with the anisotropic SU(3) lattice QCD at finite temperature below $`T_c`$, and have presented our preliminary results. We have seen the tendency that the peak becomes broader with increasing temperature, which supports the thermal width broadening  of the glueball near the critical temperature. It is interesting to apply MEM to the glueball correlator above $`T_c`$, where it becomes less trivial to find a proper ansatz for the spectral function. Note that there may survive some non-perturbative effects even above the critical temperature. In fact, it was suggested, in $`\sigma `$-$`\pi `$ sector, that a strong correlation may still survive above the critical temperature . To investigate the correlation above $`T_c`$ with MEM, it would be necessary to analyze the non-smeared correlators.

Acknowledgment

N.I. thanks Dr. S. Sasaki for useful information.


text:
# Flexibility of 𝛽-sheets: Principal-component analysis of database protein structures

## I Introduction

Most protein folds can be viewed as compact packings of a fixed set of secondary-structural elements: $`\alpha `$-helices and $`\beta `$-sheets. It can be reasoned that the formation of these elements greatly simplifies the folding free-energy landscape by reducing the number of degrees of freedom. As a first approximation, helices and sheets can be considered as rigid objects, possessing only six degrees of freedom each (three translations and three rotations). However, most helices and sheets display some amount of bending in a protein’s final fold. Understanding to what extent these elements are flexible, and which are their dominant degrees of freedom, will help to further our understanding of how proteins fold, and even how they function.

A number of studies have extracted the flexible motions of biological molecules using normal-mode and principal-component analysis (PCA). A recent PCA analysis of $`\alpha `$-helices from a structural database revealed three “soft” modes: two degenerate bend modes and a twist mode. For all but the longest helices, these three modes were sufficient to describe the deformations observed in real structures.

Arguably, a quantitative understanding of flexibility is more important for $`\beta `$-sheets than for $`\alpha `$-helices. In natural structures, sheets display a variety of highly distorted and bent shapes, e.g. $`\beta `$-barrels and twisted sheets, while helices are generally much less distorted. What are the dominant collective motions of $`\beta `$-sheets and how do they depend on the sheet’s size? Also how does sheet geometry affect the flexibility of a sheet? The “geometry” of a $`\beta `$-sheet is the amino-terminal to carboxyl-terminal orientation of the various strands making up the sheet. Most sheets fall into one of two geometries - parallel, where all the strands are oriented in the same direction, or anti-parallel, where strands alternate direction. The geometry dictates the hydrogen-bonding pattern within the sheet and hence plays a role in determining the sheet’s flexibility.

Here, we report a principal-component analysis of the flexibility of parallel and anti-parallel $`\beta `$-sheets from the Protein Data Bank. The sheets considered range in size from 3 to 6 strands, with 3 to 6 residues per strand. For both parallel and anti-parallel sheets, we find two dominant modes of flexibility: twisting about an in-plane axis that is perpendicular to the strand orientation, and bending of this same axis. The distributions of amplitudes for these two modes are independent Gaussians. Thus, the principal-component modes can be interpreted as dynamical normal modes of an elastic object. Motivated by this interpretation, we consider the scaling of mode eigenvalues (variances of amplitudes) with sheet size, and compare to predictions of a simple elastic model. For all sizes considered, parallel sheets are more rigid than anti-parallel sheets.

Recently, $`\beta `$-sheet structures have been characterized in detail by Ho and Curmi.. This database study focused on average properties of sheets, including twist, shear, and hydrogen bonding. In contrast, our PCA analysis of $`\beta `$-sheets provides a characterization of the flexibility of sheets about their average structures. Possible applications of the results reported here on sheet flexibility include improved parameterization of force-field models and inclusion of sheet elastic energies in $`\beta `$-protein design.

## II Results

### A Principal-component analysis of database $`\beta `$-sheets

The first step in analyzing the flexibility of $`\beta `$-sheets was to obtain a representative set of structures. Following the procedure described in Methods, we were able to create a database of $`3516`$ representative $`\beta `$-sheets of different geometries. In Fig. 1 we show the distribution of geometries (parallel, anti-parallel, etc.) for sheets consisting of at least 6 strands. For sheets that have more than 6 strands, we consider all the 6 stranded sub-sheets (e.g. an 8 stranded sheet would contain three 6 strand sheets). The strand directions are labeled ‘u’ for up and ‘d’ for down, as if the sheet were lying flat on the page. By convention, the first strand is always oriented in the up direction, i.e. with the amino-terminal at the bottom of the page and the carboxyl-terminal at the top of the page. Since each of the two outermost strands of a sheet can be equally well regarded as the first strand, the distribution of geometries in Fig. 1 includes every sheet twice. We find that the distribution of geometries is highly non-uniform, with the most frequent types being parallel (uuuuuu) and anti-parallel (ududud). Because other geometries occur so infrequently, we focus the rest of the analysis on parallel and anti-parallel sheets.

Sets of defect-free $`\beta `$-sheets of a given class were extracted from the 3516 representative $`\beta `$-sheets. Sheets are of the same class if they have the same size, $`S`$ strands each of length $`L`$, geometry, and pleatedness or corrugation (see Methods). To quantify flexibility, we performed a structural principal-component analysis (PCA) for each class of sheets, to identify the dominant collective fluctuations around the mean structure. To implement the PCA, we first computed the mean structure for each class of sheets via an iterative procedure. Starting with a randomly chosen sheet of the desired class, we aligned to it all other sheets of the same class. To align two sheets, we minimized the coordinate root mean square (crms) distance between their corresponding Cα atoms. A mean structure was then obtained by averaging the position of each Cα atom over all the aligned structures. This procedure was iterated, each time using the new mean structure as the basis for alignment, until the mean structure converged to within $`10^4`$ Å/residue. An example of a subset of structures aligned to the converged mean structure is shown in Fig. 2. For anti-parallel sheets we find that the average structure conforms to the sheared structure discussed in Ho and Curmi. They also showed that parallel sheets are less sheared and we find this to be the case for our average parallel sheets.

The second step in the principal-component analysis was to compute the structural covariance matrix for each class of sheet. A covariance matrix measures the correlation of the variation from the mean for each pair of coordinates. In our case, there were $`3SL`$ coordinates – 3 spatial directions for each of $`S\times L`$ Cα atoms. Consequently, the covariance matrix was a $`3SL\times 3SL`$ matrix, with elements $`i,j`$ defined as

$$C_{i,j}=\frac{1}{N1}\underset{m=1}{\overset{N}{}}(x_{mi}x_i)(x_{mj}x_j),$$
(1)
where $`N`$ is the number of sheets in the given class, $`x_{mi}`$ is the $`i^{th}`$ coordinate of the $`m^{th}`$ structure, and $`x_i`$ is the $`i^{th}`$ coordinate of the mean structure.

To complete the principal-component analysis, we computed the eigenvalues $`\{\lambda _q\}`$ and eigenvectors $`\{\stackrel{}{v}_q\}`$ of the covariance matrix for each class of sheets (available as Supplementary Material). The largest eigenvalues and corresponding eigenvectors represent the directions for which the data has the largest variance. These directions are the “soft” modes of the sheets, i.e. those collective deformations that appear with largest amplitude in the data set. Figure 3 shows the top 10 eigenvalues for anti-parallel sheets of size $`S=4`$ and $`L=5`$, as shown in Fig. 2. Each eigenvalue is given in units of Å2 and measures the variance of the distribution for a particular mode. Two dominant eigenvalues are evident in Fig. 3. The first mode is primarily a twist of the sheet about the in-plane axis perpendicular to the strand orientation (Fig. 4(a)). The second mode is primarily a bend of the sheet along the same axis (Fig. 4(b)).

For all the classes of sheets considered, we found two dominant soft modes. The eigenvalues (i.e. variances) of these modes are shown for different sheet sizes and geometries in Figs. 5 and 6. Figure 5 shows the scaling of the eigenvalues with the number of strands $`S`$ for sheets of fixed strand length. The bend-mode eigenvalues increase approximately as $`S^4`$, while the twist-mode eigenvalues increase more slowly with $`S`$ (fits to $`S^4`$ for the anti-parallel bend modes are shown). As a result of this difference in scaling, the eigenvalues for the twist and bend modes cross with increasing number of strands. For five or more strands, the eigenvalue for the bend mode becomes the larger, implying greater deformations of the sheet by bending compared to twisting. Figure 6 shows the scaling of the eigenvalues with the strand length $`L`$ for sheets with a fixed number of strands. The eigenvalues generally increase with strand length, scaling roughly as $`L`$ or $`L^2`$. The scaling behavior expected for pure bend and twist modes is discussed in the next section.

A feature that emerges from the scaling graphs is that the twist-mode and bend-mode eigenvalues are almost always larger for anti-parallel sheets than for parallel sheets. Since these two modes dominate deviations from the mean structure, this implies that total deformations are typically larger for anti-parallel sheets than for parallel sheets.

Next, we consider the actual distributions of amplitudes for the dominant twist and bend modes. The displacement of a given sheet from the mean structure, $`\delta \stackrel{}{x}=\stackrel{}{x}\stackrel{}{x}`$, can be expanded in terms of the PCA eigenvectors as $`\delta \stackrel{}{x}=_qa_q\stackrel{}{v}_q`$. The amplitude $`a_q`$ is given by the projection of the displacement vector $`\delta \stackrel{}{x}`$ onto mode $`q`$. Figures 7(a and b) show the distributions of projections onto twist and bend for the $`1454`$ anti-parallel sheets with $`S=4`$ and $`L=5`$ (cf. Fig. 2). The two distributions can be fit well by Gaussians, with the variances of the Gaussians, $`9.6985`$ Å2 for twist and $`4.3655`$ Å2 for bend, close to the exact variances given by the mode eigenvalues, $`7.7876`$ Å2 for bend and $`4.0044`$ Å2 for twist. By construction, all PCA modes are uncorrelated to lowest order, i.e. $`a_qa_q^{}=0`$, for $`qq^{}`$. To look for possible higher-order correlations, we made a scatter plot of the amplitudes for the two dominant modes, as shown in Fig. 7(c). The distributions of points is roughly ellipsoidal, indicating that there are no strong higher-order correlations between modes. Similar independent Gaussian distributions were obtained for all classes of sheets for which PCA analysis was performed.

The Gaussian distributions of mode amplitudes and the lack of higher-order correlations between modes suggest that the PCA modes can be interpreted as the dynamical modes of a $`\beta `$-sheet. This is consistent with previous results for $`\alpha `$-helices showing the near identity between PCA modes obtained from static structures and the elastic normal modes of a model helix. For small amplitude motions, elastic normal modes at equilibrium have independent Gaussian distributions $`P(a_q)`$ determined by Boltzmann weights

$$P(a_q)e^{E(a_q)/k_BT}e^{c_qa_q^2/2k_BT},$$
(2)
where $`E(a_q)=c_qa_q^2/2`$ is the deformation energy as a function of the mode amplitude $`a_q`$, and $`c_q`$ is the spring constant for the mode. The “soft” modes have the smallest spring constants, and therefore the broadest spread of amplitudes.

### B Scaling of the PCA modes

Guided by the interpretation of the PCA modes as the elastic normal modes of a sheet, we consider the scaling with sheet size of the mode eigenvalues. Let us first consider the dominant bend mode, as shown in Fig. 4(b). For a uniform bend of the in-plane axis perpendicular to the strand orientation, the displacement of the strand at position $`x`$ along this axis goes as $`\delta zx^2/R`$, where $`R`$ is the radius of curvature. The bending eigenvalue is given by $`\lambda _{\mathrm{bend}}=|\delta \stackrel{}{z}\stackrel{}{v}|^2`$ where $`\stackrel{}{v}`$ is the normalized eigenvector. It follows that

$$\lambda _{\mathrm{bend}}LS^5\frac{1}{R^2},$$
(3)
where $`L`$ is the strand length and $`S`$ is the number of strands. At thermal equilibrium, each normal mode has $`k_BT/2`$ of potential energy. For the bend mode, this energy would be put into curvature of the axis:

$$\frac{1}{2}k_BT=\frac{1}{2}\kappa LS\frac{1}{R^2},$$
(4)
where $`\kappa `$ is the bend stiffness per unit length, indicating

$$\lambda _{\mathrm{bend}}\frac{k_bT}{\kappa }S^4.$$
(5)
Thus the eigenvalue for a pure bend mode would scale as $`S^4`$ and would be independent of $`L`$. In Fig. 5, the predicted $`S^4`$ scaling of the bend-mode eigenvalue is seen for both parallel and anti-parallel $`\beta `$-sheets. In Fig. 6, however, we observe a significant increase of the eigenvalue with $`L`$, which is not predicted. The reason for the increase of the bend-mode eigenvalue with strand length $`L`$ is likely to be the significant bending of individual strands associated with this mode, which would contribute a term to $`\lambda _{\mathrm{bend}}`$ that scales as $`L^4`$.

For the twist mode, as shown in Fig. 4(a), we assume that each strand rotates by an angle $`\delta \theta `$ with respect to its neighboring strand. This corresponds to a uniform twist of the sheet about the in-plane axis perpendicular to the strand orientation. The displacement of a Cα atom on the strand at $`x`$ along this axis, and at a distance $`l`$ from this axis, is $`\delta zl(x/d)\delta \theta `$, where $`d`$ is the distance between neighboring strands. It is straightforward to show that eigenvalue for the twist mode goes as

$$\lambda _{\mathrm{twist}}\delta \theta ^2L^3S^3$$
(6)
At thermodynamic equilibrium, using $`k_BT=c\delta \theta ^2LS`$, where $`c`$ is a twist stiffness per unit length, we find $`\lambda _{\mathrm{twist}}L^2S^2`$. As shown in Fig. 6, the twist-mode eigenvalue does scale as $`L^2`$ as predicted, at least for sheets of up to $`S=5`$ strands. As seen in Fig. 5, however, the scaling of the twist-mode eigenvalue with strand number is much weaker than the predicted $`S^2`$, for all strand lengths.

## III Discussion

Our principal-component analysis of $`\beta `$-sheets indicates that sheets in proteins are deformed primarily in two ways, by twisting and bending. The amplitudes of these modes were found to have independent Gaussian distributions, suggesting that twist and bend are the “soft” elastic modes of sheets. The interpretation of the PCA twist and bend modes as elastic modes of sheets is consistent with previous PCA results on $`\alpha `$-helices. For helices, the dominant modes found by PCA were shown to be indistinguishable from the soft elastic modes of a model helix.

In this light, the generally larger twist and bend deformations found for anti-parallel sheets compared to parallel sheets suggests that anti-parallel sheets have softer elastic spring constants. This could arise from boundary effects imposed by the differences in connectivity. Anti-parallel sheets tend to be connected by short loops which allows the sheet to be easily bent. In contrast, parallel sheets have more complex connectivity (long loops) which could impinge on their ability to deform. The stiffness may also arise due to slight differences in the hydrogen bonding pattern between parallel and sheared anti-parallel sheets. It will be interesting to see if physical modeling can capture the apparent differences in elasticity between the two types of sheets.

While the PCA modes of all $`\beta `$-sheets studied had the characteristics of elastic normal modes (i.e., independent Gaussian distributions of amplitudes) significant deviations from simple scaling with sheet size were observed. These deviations can probably be attributed to two main causes. First, the actual soft modes are not the pure twist and bend modes of rigid strands as assumed in the scaling analysis. Significant deformations do occur within individual strands - the bend mode also contains bending along the strands which would contribute to the scaling of its eigenvalue with strand length $`L`$. Second, $`\beta `$-sheets, as quasi-two-dimensional objects, are subject to strong boundary effects. For example, the constraint of global connectivity of the sheets and the need to form backbone hydrogen bonds to the outer strands are likely to introduce size-dependent effects beyond the simple scaling analysis. In contrast, the nearly ideal scaling of the dominant modes of $`\alpha `$-helices probably reflects the much weaker effect of boundaries on quasi-one-dimensional objects.

Our PCA results provide a set of scaling behaviors and structural properties that can be used to test and refine energetic models of $`\beta `$-sheets. We have found that a simple spring model, similar to one that was able to capture the normal modes of helices, did not describe well the scaling properties of sheets. This suggests that the energetics governing $`\beta `$-sheet flexibility is more complex than in helices. More detailed force fields are required and these will help to further clarify the microscopic interactions governing $`\beta `$-sheet structure and flexibility.

Another possible application of the PCA results is to protein design. A recent design scheme focussed on the packing of a fixed set of structural elements to explore the space of potential novel folds. In that work, only rigid helices were considered for the structural building blocks. Because of the greater flexibility of sheets compared to helices, extending the packing scheme to sheets will require including energies of elastic deformation. The soft modes of sheets can be easily incorporated into the design of low-energy backbone structures using the elastic energies $`E(a_q)=c_qa_q^2/2`$. Since we found only two dominant modes for sheets, elasticity can be added to the packing scheme with only two extra degrees of freedom per sheet. The soft modes could also be incorporated into models that analyze how a specific protein’s structure is changed when the sequence is mutated or redesigned. Inclusion of soft modes may prove particularly useful in redesign of binding sites, which is currently limited to the rigid-backbone approximation.

In summary, this work has used database protein structures to reveal the flexible motions of $`\beta `$-sheets. The effective spring constants and eigenvectors for the sheets studied, as well as the mean structure coordinates, are available as Supplementary Material.

## IV Methods

We compiled a set of $`\beta `$-sheet structures from $`2860`$ representative protein folds in the FSSP database. The representative structures in the FSSP are the structurally distinct folds that result from doing an all-against-all structure clustering of protein folds that have less than 25% sequence identity. Thus the set is designed to minimize fold redundancy. Making use of the structural annotations in the PDB files for the representative set, we extracted the $`\beta `$-sheets from each fold. Only sheets within a certain size range were considered: sheets had to have at least 3 strands and no more than 25, with between 3 and 15 residues per strand. Using these criteria we were able to extract $`3516`$ representative $`\beta `$-sheets.

Many of the sheets in this representative set were found to have defects. A defect occurs when there is a gap in the hydrogen-bonding pattern, or the strands are of different lengths and one overhangs the other. These two types of defects are schematically illustrated in Fig. 8. We wished to eliminate sheets with defects before performing a principal-component analysis of flexibility. To systematically identify defects, we developed a procedure for finding the optimal pairing of Cα atoms on adjacent strands. A defect is indicated when a Cα atom is left unpaired. To find the optimal pairing, we computed a distance matrix $`d_{i,j}`$ for each pair of neighboring strands,

$$d_{i,j}=|\stackrel{}{r}_i\stackrel{}{r}_j|$$
(7)
where $`\stackrel{}{r}_i`$ is the position of the $`i`$th Cα atom on the first strand and $`\stackrel{}{r}_j`$ is the position of the $`j`$th Cα atom on the second strand. The two strands could be of different lengths. We defined the optimal pairing $`ij`$ between strands to be the one which minimized

$$D=6\AA N_{\mathrm{gaps}}+\underset{\mathrm{pairs}i,j}{}d_{i,j}$$
(8)
where the last term represents a penalty of +6Å for each gap. A gap occurs when a $`C_\alpha `$ atom on one strand is not paired. In practice, we found that a gap penalty of 6Å identified gaps in good agreement with a visual inspection of the sheets. To efficiently find the optimal pairing, i.e. the one which minimizes $`D`$, we employed the Needleman-Wunsch method for global sequence alignment. We applied the alignment procedure for each neighboring pair of strands in the sheet (strands in the sheet tend to be annotated in the order they occur in the sheet, hence the need to only align neighbors). This gave us the optimal alignment of all the strands in the sheet, and allowed us to identify sheet defects.

After performing the pairwise alignment of all the strands making up a sheet we then extracted sheets of a fixed size, e.g. $`S`$ strands each of length $`L`$. This is illustrated in Fig. 8(a). We scanned the aligned sheet using a window of the specified dimension. If the window spanned a region that did not contain a gap, the positions of the Cα atoms were recorded to a file. For each extracted sheet, we also recorded the sheet geometry, i.e. the relative amino-terminal to carboxyl-terminal directions of all the strands in the sheet. The direction of the first strand was defined to be up to simplify the recording of sheet geometries. The choice of which of the two outermost strands to call the first strand was arbitrary, and this was compensated for later. Another quantity recorded for each sheet was its “pleatedness”. A $`\beta `$-sheet is corrugated - it has alternating ridges and valleys. If the first row of atoms from the strands making up a sheet comprise a “ridge” we consider the sheet to have positive pleatedness. If the first row of atoms comprise a “valley” we consider the sheet to have negative pleatedness. The two types of pleatedness are shown in Fig. 8(b). We consider sheets of the same size, geometry, and pleatedness to constitute a “class”. Only sheets of the same class are aligned for subsequent principal-component analysis (see Results).

We note that each extracted ungapped sheet in fact generates two sheets because either of the two outermost strands could be called the “first”. Having arbitrarily chosen a first strand and called its direction up, there are two possible symmetry operations to obtain the second sheet. If the last strand is also up, we flip the sheet (this causes the pleatedness to change sign). If on the other hand the last strand is oriented down, we perform a clockwise rotation of the sheet by 180o. These symmetry operations effectively double our database of sheet structures.