format: text

text:
# Pairing of Lyapunov Exponents for a Hard-Sphere Gas under Shear in the Thermodynamic Limit

## I Introduction

Non-equilibrium molecular dynamics (NEMD) simulations of Navier-Stokes equations have been used to study the shear viscosity properties of fluids for a long time. To study the coefficient of shear viscosity in Navier-Stokes equations, a carefully chosen periodic boundary condition in NEMD simulations is enough to drive the system out of equilibrium. Based on these ideas, in the beginning days of the development of this subject, an algorithm was constructed from simple Newtonian equations of motion using a periodic boundary conditions, the so-called Lees-Edwards boundary conditions . However, it was soon realized that in the absence of an explicit dependence on the shear field in this algorithm, one could not make an appropriate connection with the Green-Kubo relations, and therefore, it was difficult to deal with from a theoretical point of view. As a remedy, some other algorithm with an explicit dependence on the shear field was called for, and the DOLLS and the SLLOD algorithms were born.

The main idea behind the DOLLS and the SLLOD algorithm is an explicit dependence on the shear field, $`\gamma `$. The DOLLS algorithm was implemented first . The SLLOD equations of motion were proposed soon after , and are now preferred because they are equivalent to the boundary driven method . Both algorithms have to be supplemented by a thermostat, which continuously removes the energy generated due to the work done on the system by the shear field such that an non-equilibrium steady state (NESS), homogeneous in space, can be reached.

In this paper, we will look at the SLLOD equations of motion for a gas of hard spheres from the point of view of dynamical systems. The Lyapunov exponents of a system of particles, obeying the SLLOD equations of motion and mutually interacting by means of WCA potential with an iso-kinetic Gaussian thermostat, was first studied by Morriss . The study showed that the shear viscosity can be obtained from the sum of all the Lyapunov exponents. The simulation results in Ref.  also indicated that once the Lyapunov exponents are arranged in ascending order of magnitude, the sum of the largest and the smallest, the sum of the second largest and the second smallest and so on, were the same. The phenomenon of such pairing of the Lyapunov exponents is known as the Conjugate Pairing Rule, or the CPR. Since it is in general difficult to calculate all the individual Lyapunov exponents of a system, an extensive theoretical study soon ensued to understand the CPR for the Lyapunov exponents of systems obeying the SLLOD equations of motion. Evans and co-workers investigated this point  to conclude that the Lyapunov exponents pair exactly for general inter-particle potentials and all $`\gamma `$. In a follow-up work, Sarman et al.  carried out simulation studies in support of Ref. .

In the next few years, the connection between the dynamical systems theory and statistical mechanics saw a surge of interest. Some situations were found, where it was possible to prove that the CPR is satisfied exactly . The status of the CPR for the SLLOD and the DOLLS dynamics was revisited. For a system of particles obeying the SLLOD and the DOLLS dynamics with a WCA inter-particle potential and arbitrary $`\gamma `$, CPR was reported to be violated on the basis of simulation results, but recently it was shown that this claim is based on an erroneous analysis of the Lyapunov spectrum. However, for these two systems, no attempt of a theoretical understanding about the nature of an approximate CPR has been carried out so far. In this paper, we address and attempt to clarify these issues. We find that for a dilute gas of hard spheres obeying the SLLOD dynamics, where the masses and the radii of the spheres are not necessarily the same, and the total peculiar kinetic energy is kept constant by applying the iso-kinetic Gaussian thermostat in between collisions, the CPR is violated at the most at $`O(\gamma ^4)`$, for small $`\gamma `$, in the thermodynamic limit.  Our analysis is based on the key idea that the coefficient of friction representing the iso-kinetic Gaussian thermostat for a dilute gas of particles mutually interacting by means of a short-ranged potential and obeying the SLLOD dynamics with a small shear rate, in the NESS, reaches a fixed value in the thermodynamic limit, with $`1/\sqrt{N}`$ fluctuations, where $`N`$ is the number of particles. .

The structure of the paper is the following: in Sec. II, we describe equations of motion for the SLLOD dynamics, define the Lyapunov exponents, and discuss the sufficient conditions for an exact CPR. In Sec. III, we demonstrate how the coefficient of friction representing the iso-kinetic Gaussian thermostat is expected to reach a fixed value at the thermodynamic limit, with $`1/\sqrt{N}`$ fluctuations. In Sec. IV, we present the explicit calculations and discuss the status of an approximate CPR. To make the calculations In Secs. II-IV simple, we assume that each of the gas particles has a unit mass. Finally, we end this paper with discussions on possible generalizations, including the generalization to the case when the masses of the gas particles are arbitrary, in Sec. V.

## II The SLLOD Equations of Motion from a Viewpoint of Dynamical Systems

### A Equations of Motion

The SLLOD equations of motion describe the dynamics of a collection of $`N`$ particles constituting a fluid with a macroscopic velocity field $`𝐮=\gamma y\widehat{𝐱}`$ (i.e., the gradient of the $`x`$-component of the macroscopic fluid velocity $`𝐮`$ in the $`y`$-direction is $`\gamma `$). For simplicity, each gas particle is assumed to have a unit mass. The specific form of the SLLOD equations of motion for the $`i`$-th particle, in terms of its position $`𝐫_i`$ and peculiar momentum $`𝐩_i`$, is given by

$`\dot{𝐫_i}=𝐩_i+\gamma y_i\widehat{𝐱},\dot{𝐩_i}=𝐅_i\gamma p_{iy}\widehat{𝐱}\alpha 𝐩_i.`$ (1)
where $`𝐅_i`$ is the force on the $`i`$-th particle due to the other particles in the system. The peculiar velocity of a particle is defined as its velocity with respect to the velocity of the flow at its location and the peculiar momentum of a particle is the product of its mass and its peculiar velocity. The value of $`\alpha `$, the coefficient of friction representing the iso-kinetic Gaussian thermostat in Eq. (1), is chosen such that the total peculiar kinetic energy of the system, $`_ip_i^2/2`$, is a constant of motion, i.e.,

$`\alpha ={\displaystyle \frac{_{i=1}^N(𝐅_i𝐩_i\gamma p_{ix}p_{iy})}{_{i=1}^Np_i^2}}.`$ (2)
The SLLOD equations of motion, without the dissipative term $`\alpha 𝐩_i`$, cannot be derived from a Hamiltonian (unlike the DOLLS equations of motion).

We will use the equations of motion exclusively in terms of the particles’ positions $`𝐫_i`$ and laboratory velocity $`𝐯_i`$. This introduces the change of variable from $`𝐩_i`$ to $`𝐯_i=𝐩_i+\gamma y_i\widehat{𝐱}`$ in Eq. (1), which can then be written as

$`\dot{𝐫_i}=𝐯_i,\dot{𝐯_i}=𝐅_i+\alpha \gamma y_i\widehat{𝐱}\alpha 𝐯_i.`$ (3)
In the present context, the gas particles are hard spheres of arbitrary radii. This reduces the dynamics of the gas particles to an alternating sequence of flight segments and instantaneous binary collisions. During a flight, the dynamics of the gas particles is therefore described by Eqs. (2-3) with $`𝐅_i=0`$. As far as the dynamics in collisions is concerned, it is possible to derive the limiting behavior of the iso-kinetic Gaussian thermostat as $`𝐅_i\mathrm{}`$, but this leads to rather complicated collision rules. For the purpose of simplicity, in this paper, we choose to apply the iso-kinetic thermostat only in between collisions. Thus, at an instantaneous collision between the $`i`$-th and the $`j`$-th sphere ($`i,j=1,2,..,N;ij`$), the post-collisional positions and laboratory momenta ($`+`$ subscripts) are related to their pre-collisional values ($``$ subscripts) by

$`𝐫_{i+}`$ $`=`$ $`𝐫_i,𝐫_{j+}=𝐫_j,`$ (4)
$`𝐯_{i+}`$ $`=`$ $`𝐯_i\{(𝐯_i𝐯_j)\widehat{𝐧}_{ij}\}\widehat{𝐧}_{ij}`$ (5)
and

$`𝐯_{j+}`$ $`=`$ $`𝐯_j+\{(𝐯_i𝐯_j)\widehat{𝐧}_{ij}\}\widehat{𝐧}_{ij},`$ (6)
while the positions and the velocities of the rest of the spheres remain unchanged. Here, $`\widehat{𝐧}_{ij}`$ is the unit vector along the line joining the center of the $`i`$-th sphere to the $`j`$-th sphere at the instant of collision. Note that although in any particular collision, the peculiar kinetic energy changes over a collision, these changes are random, both in magnitude and sign, due to the randomness of the collision parameters, and hence it is quite likely that the system would reach a steady state, where the average change of peculiar kinetic energy would be zero.

To study the SLLOD dynamics as a dynamical system in three dimensions (the dimensionality does not affect our analysis), we form the $`3N`$-dimensional vectors $`𝐑=(𝐫_1,𝐫_2,\mathrm{},𝐫_N)`$, $`𝐕=(𝐯_1,𝐯_2,\mathrm{},𝐯_N)`$ and $`\widehat{𝐍}_{ij}`$, whose $`l`$-th entry is given by $`\widehat{𝐍}_{ij}^l=(\delta _{l,i}\delta _{l,j})\widehat{𝐧}_{ij}/\sqrt{2}(l=1,2,..,N)`$. Using these new variables, we write the SLLOD equations of motion in the $`(𝐫_i,𝐯_i)`$ coordinates of $`N`$ hard spheres during a flight, Eq. (3), in a compact form

$`\dot{𝐑}=𝐕,\dot{𝐕}=\alpha \gamma \mathrm{𝐂𝐑}\alpha 𝐕.`$ (7)
Here, $`𝐂`$ is a $`3N\times 3N`$ matrix with $`N\times N`$ entries, each of which is a $`3\times 3`$ matrix. In terms of the entry index $`(l,m)`$, in the $`xyz`$-basis, $`𝐂_{lm}=𝐜\delta _{lm}`$ ($`l,m=1,2,..,N`$) and

$`𝐜=\widehat{𝐱}\widehat{𝐲}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right].`$ (11)
At a collision between the $`i`$-th and the $`j`$-th sphere, the equations of motion are given by

$`𝐑_+=𝐑_{},𝐕_+=𝐕_{}2(𝐕_{}\widehat{𝐍}_{ij})\widehat{𝐍}_{ij}.`$ (12)
In our analysis hereafter, except for Sec. III, we will use only Eqs. (7-12) to describe the dynamics.

At this point, we introduce the following notations. A phase space point can be denoted as $`𝚪=(𝐑,𝐕)`$. A linear transformation on phase space can be given as a $`6N\times 6N`$ matrix. Any such matrix $`𝐏`$ can be split in terms of four $`3N\times 3N`$ sub-blocks, for which we use the notation $`𝐏^{[1]}`$, $`𝐏^{[2]}`$, $`𝐏^{[3]}`$ and $`𝐏^{[4]}`$, such that

$$𝐏=\left[\begin{array}{cc}𝐏^{[1]}& 𝐏^{[2]}\\ 𝐏^{[3]}& 𝐏^{[4]}\end{array}\right].$$
(13)
Each sub-block $`𝐏^{[i]}`$ ($`i=1,\mathrm{},4`$) itself can be divided in $`3\times 3`$ sub-blocks again, where each sub-block can be identified by two indices $`l`$ and $`m`$, $`l`$ along the horizontal and $`m`$ along the vertical direction ($`l,m=1,2,..,N`$). Such a sub-block is denoted by $`𝐏_{lm}^{[i]}`$.

### B Lyapunov Exponents for Hard-sphere Systems

To calculate the Lyapunov exponents for hard-sphere systems, let us say that the system starts at time $`t_0`$ at a phase-space location $`𝚪(t_0)(𝐑(t_0),𝐕(t_0))`$. Under time evolution, $`𝚪(t)`$ follows a trajectory in the $`6N`$-dimensional phase space, which we call the “reference trajectory”. The set of $`N`$ hard spheres would suffer a sequence of binary collisions on this trajectory. We also consider an infinitesimally displaced trajectory in the phase space, which starts at the same time $`t_0`$, but at $`𝚪^{}(t_0)=𝚪(t_0)+\delta 𝚪(t_0)`$. Under time evolution, $`𝚪^{}(t)`$ follows another trajectory, always staying infinitesimally close to the reference trajectory. This trajectory we call the “adjacent trajectory”. We also assume that the set of $`N`$ hard spheres on the reference and the adjacent trajectories suffer the same sequence of binary collisions. We denote the time evolution of the infinitesimal $`6N`$-dimensional tangent vector $`\delta 𝚪(t)`$ over time $`(tt_0)`$ by the $`6N\times 6N`$ matrix $`𝐋(tt_0)`$, i.e.,

$`\delta 𝚪(t)=𝐋(tt_0)\delta 𝚪(t_0).`$ (14)
The Lyapunov exponents are the possible exponential growth rates in time of $`|𝐋\widehat{𝚪}|`$ for different directions of unit vectors $`\widehat{𝚪}`$. We have to define the norm in an appropriate way. Making the time it takes for a sphere with a typical velocity $`v_0`$ to cross the distance of a typical radius of a sphere $`a_0`$ our unit of time (i.e., $`a_0/v_0`$ is set to $`1`$) solves the problem that the components of $`\widehat{𝚪}`$ have different dimensions. For the inner product between two tangent vectors $`\delta 𝚪^{(1)}=(\delta 𝐑^{(1)},\delta 𝐕^{(1)})`$ and $`\delta 𝚪^{(2)}=(\delta 𝐑^{(2)},\delta 𝐕^{(2)})`$ we use

$`\delta 𝚪^{(1)}|\delta 𝚪^{(2)}={\displaystyle \underset{i=1}{\overset{N}{}}}(\delta 𝐫_i^{(1)}\delta 𝐫_i^{(2)}+\delta 𝐯_i^{(1)}\delta 𝐯_i^{(2)}).`$

The norm is now defined as $`|\delta 𝚪|=\sqrt{\delta 𝚪|\delta 𝚪}`$.

The Lyapunov exponents are the logarithms of the eigenvalues of the matrix $`𝚲`$, defined by

$`𝚲=\underset{t\mathrm{}}{lim}\left[\stackrel{~}{}L(tt_0)\right]^{1/[2(tt_0)]},`$ (15)
where, $`\stackrel{~}{}L(tt_0)=[𝐋(tt_0)]^\text{T}𝐋(tt_0)`$. The corresponding directions for the exponential expansion and contraction of the phase space $`(𝐑,𝐕)`$ are obtained from the eigenvectors of $`\stackrel{~}{}L(tt_0)`$.

The dynamics of $`\delta 𝚪(t)`$ in Eq. (14), for a gas of hard spheres, can be decomposed into an alternating sequence of flights and instantaneous binary collisions. We denote the transformation of $`\delta 𝚪(t)`$ over a flight segment between $`t`$ and $`t+\mathrm{\Delta }t`$ by $`𝐇(\mathrm{\Delta }t)`$ such that

$`\delta 𝚪(t+\mathrm{\Delta }t)=𝐇(\mathrm{\Delta }t)\delta 𝚪(t)\text{with}𝐇(0)=𝐈.`$ (16)
Explicitly, $`𝐇(\mathrm{\Delta }t)`$ is obtained from

$`\dot{\delta 𝚪}(t)=𝐓(t)\delta 𝚪(t)`$ (17)
as

$`𝐇(\mathrm{\Delta }t)=\mathrm{exp}_{\text{ T}}\left[{\displaystyle _t^{t+\mathrm{\Delta }t}}𝑑t^{}𝐓(t^{})\right],`$ (18)
where the subscript T indicates time ordering. Notice that $`𝐇(\mathrm{\Delta }t)`$ in a general system will depend on time $`t`$ as well, but we have suppressed that in our notation. If we now denote the transformation of $`\delta 𝚪(t)`$ over an instantaneous binary collision (say, between the $`i`$-th and the $`j`$-th sphere) by the matrix $`𝐌_{ij}`$, we can express the matrix $`𝐋(tt_0)`$ in terms of the $`𝐇`$ and $`𝐌_{ij}`$ matrices in the following way: if the dynamics involves flight segments separated by $`s`$ instantaneous binary collisions at $`t_1`$, $`t_2\mathrm{}t_s`$ such that $`t_0<t_1<t_2<\mathrm{}<t_s<t`$, then

$`𝐋(tt_0)=𝐇(\mathrm{\Delta }t_s)𝐌_{i_sj_s}𝐇(\mathrm{\Delta }t_{s1})\mathrm{}`$ (20)
$`\mathrm{}𝐌_{i_1j_1}𝐇(\mathrm{\Delta }t_0).`$
Here, $`\mathrm{\Delta }t_i=t_{i+1}t_i`$ for $`i=1,2,\mathrm{},(s1)`$ and $`\mathrm{\Delta }t_s=tt_s`$.

### C The Sufficient Conditions for an Exact CPR

If the CPR is exactly satisfied for a dynamical system, the sum of the conjugate pairs of the Lyapunov exponents is some constant $`c`$, i.e., if $`\lambda _i`$ is a Lyapunov exponent of this system, then $`c\lambda _i`$ is also a Lyapunov exponent. The proof of a possible conjugate pairing rule will follow from the properties of the matrix $`𝐋(tt_0)`$. However, to understand the interplay between the properties of the matrix $`𝐋(tt_0)`$ and an exact CPR in full generality, below we first look at the property of $`𝐋(tt_0)`$ that has been used in various cases to prove CPR.

1. If the matrix $`𝐋(tt_0)`$ is symplectic, i.e., $`𝐋(tt_0)`$ satisfies the symplectic condition

$`[𝐋(tt_0)]^\text{T}\mathrm{𝐉𝐋}(tt_0)=𝐉,`$

with $`𝐉`$ as the usual symplectic matrix, then

$`\stackrel{~}{}L(tt_0)𝐉\stackrel{~}{}L(tt_0)=𝐉.`$ (21)
Eq. (21) can be used to show that

$`\text{Det}[\stackrel{~}{}L(tt_0)\stackrel{~}{L}𝐈]=0`$ (22)
(with $`𝐈`$ the identity matrix) implies

$`\text{Det}[\stackrel{~}{}L(tt_0){\displaystyle \frac{1}{\stackrel{~}{L}}}𝐈]=0.`$ (23)
This means that if $`\stackrel{~}{L}`$ is an eigenvalue of $`\stackrel{~}{}L(tt_0)`$ then so is $`\stackrel{~}{L}^1`$. Since the Lyapunov exponents are the logarithms of the eigenvalues of the matrix $`𝚲`$, defined in Eq. (15), it is easy to see that $`c=0`$. All Hamiltonian systems fall in this class.
2. In the existing literature , the concept of the symplectic condition defined above, has been generalized to the so-called “$`\mu `$-symplectic condition” and applied to thermostatted systems where an iso-kinetic Gaussian thermostat keeps the total laboratory kinetic energy constant and the external force on the constituent particles of the system is dependent only on the positions of the particles. For these systems, in an appropriate reduced phase space characterized by all the non-zero Lyapunov exponents, the matrix $`𝐋(tt_0)`$ satisfies this $`\mu `$-symplectic condition, which means that there exists a time-dependent positive scalar quantity $`\mu `$, such that $`[𝐋(tt_0)]^\text{T}\mathrm{𝐉𝐋}(tt_0)=\mu 𝐉`$. This implies that

$`\stackrel{~}{}L(tt_0)𝐉\stackrel{~}{}L(tt_0)=\mu ^2𝐉;`$ (24)
which can be used to derive that if

$`\text{Det}[\stackrel{~}{}L(tt_0)\stackrel{~}{L}𝐈]=0`$ (25)
then also

$`\text{Det}[\stackrel{~}{}L(tt_0){\displaystyle \frac{\mu ^2}{\stackrel{~}{L}}}𝐈]=0,`$ (26)
for an eigenvalue $`\stackrel{~}{L}`$ of the matrix $`\stackrel{~}{}L(tt_0)`$. That means that if $`\stackrel{~}{L}`$ is an eigenvalue of $`\stackrel{~}{}L(tt_0)`$ then so is $`\mu ^2\stackrel{~}{L}^1`$. In that case, one finds from Eq. (15) that $`c=lim_t\mathrm{}(\mathrm{ln}\mu )/(tt_0)`$. If the system is ergodic, then this long time average for $`c`$ can be equated to a NESS average. Notice that condition (a) is obtained as special case of condition (b), namely when $`\mu =1`$.

Returning momentarily to the SLLOD dynamics, we observe that the formalism developed in Refs.  fails here. The primary reason is associated with the fact that the total peculiar kinetic energy is held constant for the SLLOD dynamics, as opposed to the total laboratory kinetic energy in Refs. . One however needs to interpret the statement regarding the connection between the violation of condition (b) and the non-exactness of the CPR with care. By virtue of the fact that condition (b) above is a sufficient condition for the CPR to hold exactly, the violation of an exact CPR cannot be guaranteed if condition (b) is not satisfied.

Guided by this observation, the interplay between the properties of the matrix $`𝐋(tt_0)`$ and an exact CPR for a dynamical system can be generalized further than what is presented in (b). If there exists any constant non-singular matrix $`𝐊`$ satisfying $`𝐊^2𝐈`$ and the following condition

$`[𝐋(tt_0)]^\text{T}𝐊𝐋(tt_0)=\mu 𝐊`$ (27)
is satisfied with a time-dependent scalar quantity $`\mu `$, then Eq. (26) can be shown to hold for an eigenvalue $`\stackrel{~}{L}`$ of the matrix $`\stackrel{~}{}L(tt_0)`$, implying that the CPR is exactly satisfied for such a dynamical system . In analogy with the nomenclature presented in (b), we call Eq. (27) a “generalized $`\mu `$-symplectic condition” with matrix $`𝐊`$. We emphasize that the necessary condition for an exact CPR to hold for a dynamical system is not known.

In view of Eq. (27), thus, one should look for such a matrix $`𝐊`$ to prove an exact CPR. Instead, we look at it from a different angle, namely that we would like to understand how the SLLOD dynamics of $`N`$ hard spheres with an iso-kinetic Gaussian thermostat deviates from an exact CPR.

## III Behavior of $`\alpha `$ in the Thermodynamic Limit

Our procedure to study the deviations from CPR begins with the following observation: in the thermodynamic limit, for the SLLOD dynamics with short-range inter-particle potentials at low density of spheres and at small $`\gamma `$, the behavior of $`\alpha `$ simplifies to a great extent. After some transient time the system reaches the NESS, and the coefficient of friction $`\alpha `$ fluctuates with $`1/\sqrt{N}`$ fluctuations around a fixed value $`\alpha _0`$ . For not too large fluctuations, the distribution function for $`\alpha `$ can also be shown to be approximately Gaussian. Thus, to calculate the Lyapunov exponents for large $`N`$ at low density of spheres and at small $`\gamma `$, to which we confine ourselves henceforth, $`\alpha `$ can be replaced by $`\alpha _0`$ in Eq. (7). We will now briefly present the gist of the derivation in Ref. , applied to hard spheres.

For a hard sphere system, the force term in Eqs. (1) and (2) is zero during a flight. Thus, for a flight, we have

$`\alpha =\gamma \left[{\displaystyle \underset{i=1}{\overset{N}{}}}p_i^2\right]^1{\displaystyle \underset{i=1}{\overset{N}{}}}p_{ix}p_{iy}.`$ (28)
Introducing a second thermostat variable,

$`\beta =\gamma ^2\left[{\displaystyle \underset{i=1}{\overset{N}{}}}p_i^2\right]^1{\displaystyle \underset{i=1}{\overset{N}{}}}p_{iy}^2,`$ (29)
a closed set of equation follows from Eq. (1):

$`\dot{\alpha }`$ $`=`$ $`\mathrm{\hspace{0.17em}2}\alpha ^2+\beta `$ (30)
$`\dot{\beta }`$ $`=`$ $`\mathrm{\hspace{0.17em}2}\alpha \beta .`$ (31)
These equations are valid during the flights, i.e., the intervals between collisions.

We treat collisions by looking at their net effect, i.e., how the velocities and positions of the particles $`i`$ and $`j`$ involved in the collision, are changed from their pre-collisional values $`𝐩_i`$ and $`𝐩_j`$ to their post-collisional values $`𝐩_{i+}`$ and $`𝐩_{j+}`$. These are the only two velocities to change, and because $`_{i=1}^Np_i^2`$ is of order $`N`$, the changes in $`\alpha `$ and $`\beta `$ are, according to Eqs. (29) and (31), of order $`N^1`$. The number of collisions in the whole system is an extensive quantity as well, so there are $`O(N)`$ of these $`O(1/N)`$ changes in a unit of time. The averages of the small changes are not zero, so there is a net effect of $`O(1)`$ per unit time to the time derivatives of $`\alpha `$ and $`\beta `$, which we will denote by $`a`$ respectively, $`b`$:

$`\dot{\alpha }`$ $`=`$ $`\mathrm{\hspace{0.17em}2}\alpha ^2+\beta +a`$ (32)
$`\dot{\beta }`$ $`=`$ $`\mathrm{\hspace{0.17em}2}\alpha \beta +b.`$ (33)
This set of equations has a fixed point $`(\alpha _0,\beta _0)`$, which is stable if $`\alpha _0>0`$, so the system reaches this fixed point after some time. On top of this dynamics, there are fluctuations. Assuming that the collisions are independent, the central limit theorem applies, and the fluctuations are $`O(1/\sqrt{N})`$. For more detailed analysis, we refer to Ref. , where the independence of the changes in $`\alpha `$ and $`\beta `$ is linked to the assumption of molecular chaos.

Finally, we note that to maintain a stationary total peculiar kinetic energy in a system with a constant $`\alpha _0`$ thermostat, this constant has to be chosen differently for varying $`\gamma `$. For $`\gamma =0`$, i.e., in equilibrium, $`\alpha _0`$ can be set to zero, and kinetic energy is determined by the initial conditions. Near equilibrium, i.e., in the linear response regime, the right hand side of Eq. (28) should scale as $`\gamma ^2`$, i.e. $`\alpha _0\gamma ^2`$. Obviously, there are higher order corrections to this behavior which play a role for larger values of $`\gamma `$. If, for $`\gamma 0`$, the initial condition is such that the total peculiar kinetic energy is not equal to the stationary value, that value will be approached in time.

## IV Status of an Approximate CPR in the Thermodynamic Limit

Based on the discussion in the last paragraph of Sec. II and using the results in Sec. III, we will explore the possibility of an approximate CPR for the the SLLOD dynamics of hard spheres in the thermodynamic limit, at small $`\gamma `$ and at low density in this section. We will first obtain the desired results for a constant coefficient of friction $`\alpha _0`$ in the equations of motion (7) and then we will discuss the validity of an approximate CPR when the system is under an iso-kinetic Gaussian thermostat. To this end, our starting aim is to study the generalized $`\mu `$-symplectic properties of the matrix $`𝐋(tt_0)`$ for the dynamics described by Eq. (7) during a flight and Eq. (12) during an instantaneous binary collision between the $`i`$-th and the $`j`$-th sphere.

However, as the matrix $`𝐋(tt_0)`$ is constructed from the $`𝐇`$ and the $`𝐌`$ matrices, we will have to study the generalized $`\mu `$-symplectic properties of the $`𝐇`$ and the $`𝐌`$ matrices separately.

### A Generalized $`\mu `$-symplecticity Property of $`𝐇(\mathrm{\Delta }t)`$

The matrix $`𝐓`$ describing the dynamics of $`\delta 𝚪`$ during flights is found from Eq. (7) to be

$`𝐓\left[\begin{array}{cc}\mathrm{𝟎}& 𝐈\\ \alpha _0\gamma 𝐂& \alpha _0𝐈\end{array}\right].`$ (36)
From Eq. (36), it is straightforward to obtain

$`𝐇(\mathrm{\Delta }t)=\mathrm{exp}[𝐓\mathrm{\Delta }t],`$ (37)
where the $`3\times 3`$ sub-blocks of $`𝐇`$ are given by $`𝐇_{lm}^{[k]}(\mathrm{\Delta }t)=𝐡^{(k)}(\mathrm{\Delta }t)\delta _{l,m}`$, where

$`𝐡^{(1)}(\mathrm{\Delta }t)`$ $`=`$ $`𝐈+\left[\gamma \mathrm{\Delta }t{\displaystyle \frac{\gamma [1\mathrm{exp}(\alpha _0\mathrm{\Delta }t)]}{\alpha _0}}\right]𝐜,`$ (38)
$`𝐡^{(2)}(\mathrm{\Delta }t)`$ $`=`$ $`{\displaystyle \frac{1e^{\alpha _0\mathrm{\Delta }t}}{\alpha _0}}𝐈`$ (40)
$`+{\displaystyle \frac{\gamma }{\alpha _0^2}}\left[\alpha _0\mathrm{\Delta }t(1+e^{\alpha _0\mathrm{\Delta }t})2+2e^{\alpha _0\mathrm{\Delta }t}\right]𝐜,`$
$`𝐡^{(3)}(\mathrm{\Delta }t)`$ $`=`$ $`\gamma [1e^{\alpha _0\mathrm{\Delta }t}]𝐜\text{and}`$ (41)
$`𝐡^{(4)}(\mathrm{\Delta }t)`$ $`=`$ $`e^{\alpha _0\mathrm{\Delta }t}\{𝐈\gamma \left[\mathrm{\Delta }t+{\displaystyle \frac{1}{\alpha _0}}[1e^{\alpha _0\mathrm{\Delta }t}]\right]𝐜\}.`$ (42)
Due to the complicated form of $`𝐇(\mathrm{\Delta }t)`$, it is easier to study its $`\mu `$-symplecticity properties in terms of the matrix $`𝐓`$. This involves the task of finding a possible matrix $`𝐊`$ satisfying the condition

$`𝐓^\text{T}𝐊+\mathrm{𝐊𝐓}=\beta 𝐊`$ (43)
such that $`𝐊^2𝐈`$. If such a matrix $`𝐊`$ exists, then $`𝐇(\mathrm{\Delta }t)`$ is generalized $`\mu `$ symplectic with that matrix, and

$`\mu =\mathrm{exp}\left[{\displaystyle _t^{t+\mathrm{\Delta }t}}\beta 𝑑t^{}\right].`$ (44)
Since $`𝐊`$ and $`𝐓`$ are constant matrices in the present context, $`\beta `$ is also a constant. Eq. (43) can be treated as a simple eigenvalue equation to solve for the eigenvalue $`\beta `$ and the eigenvector $`𝐊`$. We find that there exists a matrix $`𝐆`$ satisfying

$`𝐓^\text{T}𝐆+\mathrm{𝐆𝐓}=\alpha _0𝐆`$ (45)
and the $`6N\times 6N`$ matrix $`𝐆`$ in terms of its $`3\times 3`$ sub-blocks, is given by $`𝐆^{[1]}=𝐆^{[4]}=\mathrm{𝟎}`$ and $`𝐆_{lm}^{[2]}=𝐆_{lm}^{[3]}=𝐠`$, where

$`𝐠=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right].`$ (49)
We note that there may exist other forms of $`𝐠`$ such that $`𝐆`$ satisfies Eq. (45), but Eq. (49) is the simplest one that satisfies $`𝐆^2𝐈`$ and works for all $`\gamma `$.

For the purpose of future use, we construct matrices $`𝐓_0`$ and $`𝐇_0(\mathrm{\Delta }t)`$ by setting $`\gamma =0`$ in Eq. (36) and Eq. (42), without setting $`\alpha _0=0`$ (in reality, $`\alpha _0=0`$ when $`\gamma =0`$), i.e.,

$`𝐓_0`$ $``$ $`\left[\begin{array}{cc}\mathrm{𝟎}& 𝐈\\ \mathrm{𝟎}& \alpha _0𝐈\end{array}\right]\text{and}`$ (52)
$`𝐇_0(\mathrm{\Delta }t)`$ $``$ $`\mathrm{exp}\left[𝐓_0\mathrm{\Delta }t\right].`$ (53)
More explicitly, the $`3\times 3`$ sub-blocks of $`𝐇_0(\mathrm{\Delta }t)`$ are given by $`(𝐇_0)_{lm}^{[k]}(\mathrm{\Delta }t)=𝐡_0^{(k)}(\mathrm{\Delta }t)\delta _{l,m}`$, and $`𝐡_0^{(k)}(\mathrm{\Delta }t)`$ can be found by putting $`\gamma =0`$ in Eq. (42) without putting $`\alpha _0=0`$. The matrix $`𝐇_0(\mathrm{\Delta }t)`$ is now not only generalized $`\mu `$-symplectic with matrix $`𝐆`$, but also $`\mu `$-symplectic with $`𝐉`$, i.e.,

$`[𝐇_0(\mathrm{\Delta }t)]^\text{T}𝐉𝐇_0(\mathrm{\Delta }t)=e^{\alpha _0\mathrm{\Delta }t}𝐉.`$ (54)
The relevance of this observation will become clear in Sec. IV.C.

We also note that Eqs. (36-49) hold for any constant coefficient of friction (not necessarily $`\alpha _0`$), which implies, using Eq. (27), that for a collisionless gas of point particles obeying the SLLOD dynamics with a constant coefficient of friction, the CPR is exact, as can be seen in simulations .

### B Generalized $`\mu `$-symplecticity Property of $`𝐌_{ij}`$

Unlike the $`𝐇`$ matrices, the $`𝐌_{ij}`$ matrices corresponding to a binary collision between the $`i`$-th and the $`j`$-th sphere do not follow from Eq. (12) directly. This is due to the fact that even though the sequence of binary collisions are the same on the reference and the adjacent trajectories, the binary collision between the $`i`$-th and the $`j`$-th sphere on these two trajectories in the phase space are not simultaneous. One therefore needs the dynamics of the tangent vectors for the time interval, $`\delta \tau `$, between the two collisions on the reference and the adjacent trajectories involving the $`i`$-th and the $`j`$-th sphere. To obtain an expression of $`𝐌_{ij}`$, we follow the formalism developed in Refs. , which in turn, is based on the formalism presented in Refs.  and .

The dynamics of tangent vectors at a collision is derived in the Appendix. The result is that

$`\delta 𝚪_+=𝐌_{ij}\delta 𝚪_{}`$ (55)
with

$`𝐌_{ij}=(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})\left[\begin{array}{cc}𝐈& \mathrm{𝟎}\\ 𝐑& 𝐈\end{array}\right],`$ (58)
and $`𝐑`$ a symmetric matrix given by Eq. (A29). This form of $`𝐌_{ij}`$ immediately implies that $`𝐌`$ is symplectic, but not generalized $`\mu `$-symplectic with matrix $`𝐆`$ for $`\mu =1`$, i.e.,

$`𝐌_{ij}^\text{T}\mathrm{𝐉𝐌}_{ij}=𝐉`$ (59)
but

$`𝐌_{ij}^\text{T}\mathrm{𝐆𝐌}_{ij}𝐆.`$ (60)
### C Generalized $`\mu `$-symplecticity Property of $`𝐋(tt_0)`$ and the Origin of an Approximate CPR

From Secs. IV.A and IV.B above, we can finally see that for a collection of hard spheres obeying the SLLOD equations of motion with constant coefficient of friction $`\alpha _0`$

1. the $`𝐇`$ matrices are generalized $`\mu `$-symplectic with matrix $`𝐆`$, but not with matrix $`𝐉`$ (see Eq. (45))
2. the $`𝐌`$ matrices are symplectic but not generalized $`\mu `$-symplectic with matrix $`𝐆`$ (see Eqs. (59) and (60)) .

Once the $`𝐇`$ and the $`𝐌`$ matrices are combined together, following Eq. (20), the matrix

$`𝐋(tt_0)=𝐇(\mathrm{\Delta }t_s)𝐌_{i_sj_s}𝐇(\mathrm{\Delta }t_{s1})\mathrm{}𝐌_{i_1j_1}𝐇(\mathrm{\Delta }t_0)`$ (61)
is seen to be generalized $`\mu `$-symplectic neither  with matrix $`𝐆`$ nor  with matrix $`𝐉`$.

To study the degree of deviation from an exact CPR using the properties of $`𝐋(tt_0)`$ in Eq. (61) with constant coefficient of friction $`\alpha _0`$, we can use either $`𝐊=𝐆`$ or $`𝐊=𝐉`$. While the former choice implies that one has to try to estimate the deviation from an exact CPR using the distribution of the unit vectors $`\widehat{𝐍}_{ij}`$’s and the collision angles for different sets of binary collisions in the expression of $`𝐌`$, the latter choice means that one can make the estimate by using the property of the $`𝐇`$ matrices in Eq. (61). We choose the latter approach, because not only is it much easier to calculate the typical magnitude of a flight time of a sphere at low densities, but also, an estimate of the deviation from the exact CPR can be made at small $`\gamma `$, as a power series expansion in $`\gamma `$. However, the smallness of $`\gamma `$, which has a dimension of inverse time, has to be defined in a proper manner. To do so, we notice that the density of the spheres $`n`$ sets a time scale in the form of the mean flight time $`\tau _0`$ of an individual sphere, and in three dimensions $`\tau _0\stackrel{~}{n}a_0/v_0`$. Here, $`\stackrel{~}{n}=na_0^3`$ is the dimensionless density and, as before, $`v_0`$ and $`a_0`$ are typical velocity and radius of a sphere. Thus, the actual dimensionless small parameter corresponding to the shear rate is $`\stackrel{~}{\gamma }=\gamma \tau _0`$.

A naive way to estimate the deviation from an exact CPR using the latter approach is to use the deviation of the $`𝐇(\mathrm{\Delta }t)`$ matrices from an exact $`\mu `$-symplecticity (see Eq. (61)). Such a deviation is characterized by the matrix $`𝐃(\mathrm{\Delta }t)=[𝐇(\mathrm{\Delta }t)]^\text{T}\mathrm{𝐉𝐇}(\mathrm{\Delta }t)e^{\alpha _0\mathrm{\Delta }t}𝐉`$. The matrix $`𝐃(\mathrm{\Delta }t)`$ can easily be calculated from Eq. (42). However, to estimate the order of the matrix elements of $`𝐃(\mathrm{\Delta }t)`$, an order of estimate of the quantity $`\mathrm{\Delta }t`$ has to be obtained. To this end, we note that while $`\tau _0`$ is the mean flight time for an individual sphere, $`\mathrm{\Delta }t`$ in Eq. (61) denotes the mean time for a flight of $`N`$ spheres. This implies that $`\mathrm{\Delta }t\tau _0/N`$, as on an average, there are $`N/2`$ different binary collisions over a mean flight time $`\tau _0`$ of an individual sphere. Thus, one would expect that in the thermodynamic limit, $`𝐃(\mathrm{\Delta }t)𝐃(0)=\mathrm{𝟎}`$ and one would be led to conclude that the $`𝐇`$ matrices in Eq. (61) are all symplectic. This in turn would imply, from Eq. (61), that $`𝐋(tt_0)`$ would be $`\mu `$-symplectic and therefore a gas of hard spheres, obeying the SLLOD dynamics with a constant coefficient of friction $`\alpha _0`$ would satisfy an exact CPR in the thermodynamic limit. We demonstrate below that this simplification is not correct.

The proper estimate of the deviation from an exact CPR has to be made by considering $`𝐇(\tau _0)`$. To see why this is so, we rewrite the matrix $`𝐇(\mathrm{\Delta }t)`$ as

$`𝐇(\mathrm{\Delta }t)=e^{\frac{N1}{2}\alpha _0\mathrm{\Delta }t}{\displaystyle \underset{i=1}{\overset{N}{}}}𝐇^i(\mathrm{\Delta }t)`$ (62)
with $`𝐇^i(\mathrm{\Delta }t)`$ defined by

$`(𝐇^i)_{lm}^{[k]}(\mathrm{\Delta }t)`$ $`=`$ $`\delta _{l,m}[\delta _{l,i}𝐡^{(k)}(\mathrm{\Delta }t)`$ (64)
$`+(1\delta _{l,i})(\delta _{k,1}+\delta _{k,4})e^{\alpha _0\mathrm{\Delta }t/2}𝐈].`$
In effect, $`𝐇^i(\mathrm{\Delta }t)`$ describes the evolution of the infinitesimal deviation of the trajectory of the $`i`$-th sphere, while it has an almost trivial action on the infinitesimal deviation of the trajectory of the $`j`$-th sphere, $`ji`$. It is easy to see that $`𝐇^i(\mathrm{\Delta }t)`$ has some useful properties

$`𝐇^i(\mathrm{\Delta }t_1)𝐇^i(\mathrm{\Delta }t_2)`$ $`=`$ $`𝐇^i(\mathrm{\Delta }t_1+\mathrm{\Delta }t_2),`$ (65)
$`[𝐇^i(\mathrm{\Delta }t_i),𝐇^j(\mathrm{\Delta }t_j)]`$ $`=`$ $`0\text{and}`$ (66)
$`[𝐇^i(\mathrm{\Delta }t),𝐌_{i_cj_c}]`$ $`=`$ $`0\text{if }i_ci\text{ and }j_ci.`$ (67)
The properties of $`𝐇^i(\mathrm{\Delta }t)`$ in Eq. (67) allow us to shuffle the terms in Eq. (61) such as to collect together as many $`𝐇^i(\mathrm{\Delta }t)`$s with the same $`i`$ as possible. The result is that to the right of any $`𝐌_{i_cj_c}`$ figure an $`𝐇^{i_c}(\tau _{i_c}^c)`$ and an $`𝐇^{j_c}(\tau _{j_c}^c)`$, where $`\tau _{i_c}^c`$ and $`\tau _{j_c}^c`$ are the time of flights for the $`i_c`$-th and the $`j_c`$-th spheres before their mutual collision $`c`$. Consequently,

$`𝐋(tt_0)`$ $`=`$ $`e^{\frac{N1}{2}\alpha _0(tt_0)}\left\{{\displaystyle \underset{i=1}{\overset{N}{}}}𝐇^i(t\stackrel{~}{t}_i)\right\}`$ (69)
$`\times {\displaystyle \underset{c=1}{\overset{s}{}}}𝐌_{i_cj_c}𝐇^{i_c}(\tau _{i_c}^c)𝐇^{j_c}(\tau _{j_c}^c).`$
The product sign in Eq. (69) is to be expanded towards the left, i.e. $`_{c=1}^sA_c=A_s\mathrm{}A_1`$. Here, $`\stackrel{~}{t}_i`$ is the last time that $`i`$-th particle collided (or $`t_0`$ if it didn’t collide). From Eq. (69), it is now clear that the proper estimate for the deviation from an exact CPR has to be made by considering the properties of $`𝐇^i(\tau _i)`$, with $`\tau _i=O(\tau _0)`$, and not from the properties of $`𝐇^i(\tau _0/N)`$.

We also notice that if one uses the corresponding $`𝐇_0`$ matrices instead of the $`𝐇`$ matrices in Eq. (61) to construct an analogous matrix $`𝐋_0(tt_0)`$, defined by

$`𝐋_0(tt_0)`$ $`=`$ $`𝐇_0(\mathrm{\Delta }t_s)𝐌_{i_sj_s}𝐇_0(\mathrm{\Delta }t_{s1})\mathrm{}`$ (71)
$`\mathrm{}𝐌_{i_1j_1}𝐇_0(\mathrm{\Delta }t_0),`$
then the matrix $`𝐋_0(tt_0)`$ is $`\mu `$-symplectic, because of Eqs. (54) and (59). As a consequence, the logarithms of the eigenvalues of $`\stackrel{~}{}L_0(tt_0)`$, defined by

$`\stackrel{~}{}L_0(tt_0)=[𝐋_0(tt_0)]^\text{T}𝐋_0(tt_0),`$ (72)
pair exactly. If we arrange the corresponding Lyapunov spectrum

$`𝚲_0`$ $`=`$ $`\underset{t\mathrm{}}{lim}\left[\stackrel{~}{}L_0(tt_0)\right]^{1/2(tt_0)},`$ (73)
in the decreasing order of magnitude as $`\lambda _1^{(0)}`$, $`\lambda _2^{(0)},`$ $`\mathrm{},\lambda _{6N}^{(0)}`$, then $`\lambda _i^{(0)}+\lambda _{6Ni+1}^{(0)}=\alpha _0`$.

Motivated by this, our approach to study the deviation from an exact CPR for the matrix $`𝐋(tt_0)`$ will be to take $`𝐋_0(tt_0)`$ as the reference matrix. This is reasonable because, as we will show, $`𝐋(tt_0)`$ and $`𝐋_0(tt_0)`$ are very close for small $`\stackrel{~}{\gamma }`$, if $`tt_0=\tau =O(\tau _0)`$. To show this, we first write the matrix $`𝐋_0(tt_0)`$ in the same form as Eq. (69), with $`𝐇^i(\tau _i)`$s replaced by $`𝐇_0^i(\tau _i)`$s. We then relate, for $`tt_0=\tau =O(\tau _0)`$, the difference matrix $`𝚫L(tt_0)=𝐋(tt_0)𝐋_0(tt_0)`$ to the difference between $`𝐇^i(\tau _i)`$ and $`𝐇_0^i(\tau _i)`$.

We now define new matrices $`\mathrm{\Delta }𝐇^i`$ and $`𝐄^i`$, such that

$`\mathrm{\Delta }𝐇^i(\tau _i)`$ $``$ $`𝐇^i(\tau _i)𝐇_0^i(\tau _i)`$ (74)
$`𝐄^i(\tau _i)`$ $``$ $`\left[𝐇_0^i(\tau _i)\right]^1\mathrm{\Delta }𝐇^i(\tau _i)`$ (75)
with

$`(𝐇_0^i)_{lm}^{[k]}(\tau _i)`$ $`=`$ $`\delta _{l,m}[\delta _{l,i}𝐡_0^{(k)}(\tau _i)`$ (77)
$`+(1\delta _{l,i})(\delta _{k,1}+\delta _{k,4})e^{\alpha _0\tau _i/2}𝐈]`$
\[for the definition of $`𝐡_0^{(k)}(\tau _i)`$, see the preceding paragraph of Eq. (54)\]. It is easily seen that $`𝐄^i(\tau _i)`$ has non-zero elements only for entries involving the $`i`$-th sphere. This, together with relations (67), implies that

$`[𝐄^i(\tau _i),𝐇_0^j(\tau _j)]`$ $`=`$ $`0\text{for }ij\text{ and}`$ (78)
$`[𝐄^i(\tau _i),𝐌_{i_cj_c}]`$ $`=`$ $`0\text{if }i_ci\text{ and }j_ci.`$ (79)
The matrix $`𝐄^i(\tau _i)`$ can be easily calculated from Eqs. (75). With the aid of Eqs. (42) and (53), we now express the matrix $`\mathrm{\Delta }𝐇^i(\tau _i)`$ defined in Eq. (75) as $`\mathrm{\Delta }𝐇_{lm}^{i(k)}(\tau _i)=\mathrm{\Delta }𝐡^{i[k]}(\tau _i)\delta _{l,m}`$, where

$`\mathrm{\Delta }𝐡^{i(1)}(\tau _i)`$ $`=`$ $`\left[\gamma \tau _i{\displaystyle \frac{\gamma [1\mathrm{exp}(\alpha _0\tau _i)]}{\alpha _0}}\right]𝐜,`$ (80)
$`\mathrm{\Delta }𝐡^{i(2)}(\tau _i)`$ $`=`$ $`{\displaystyle \frac{\gamma }{\alpha _0^2}}\left[\alpha _0\tau _i(1+e^{\alpha _0\tau _i})2+2e^{\alpha _0\tau _i}\right]𝐜,`$ (81)
$`\mathrm{\Delta }𝐡^{i(3)}(\tau _i)`$ $`=`$ $`\gamma [1e^{\alpha _0\tau _i}]𝐜\text{and}`$ (82)
$`\mathrm{\Delta }𝐡^{i(4)}(\tau _i)`$ $`=`$ $`\gamma e^{\alpha _0\tau _i}\left[\tau _i+{\displaystyle \frac{1}{\alpha _0}}[1e^{\alpha _0\tau _i}]\right]𝐜.`$ (83)
The crucial point is that since $`\alpha _0\gamma \stackrel{~}{\gamma }`$, i.e., the dissipation is quadratic in the shear field, $`\mathrm{\Delta }𝐇^i(\tau _i)`$ in Eq. (83) can be expanded in powers of $`\gamma `$ to obtain

$`\mathrm{\Delta }𝐡^{i(1)}(\tau _i)`$ $`=`$ $`{\displaystyle \frac{\gamma \alpha _0\tau _i^2}{2!}}𝐜,`$ (84)
$`\mathrm{\Delta }𝐡^{i(2)}(\tau _i)`$ $`=`$ $`{\displaystyle \frac{\gamma \alpha _0\tau _i^3}{3!}}𝐜,`$ (85)
$`\mathrm{\Delta }𝐡^{i(3)}(\tau _i)`$ $`=`$ $`\gamma \alpha _0\tau _i𝐜\text{and}`$ (86)
$`\mathrm{\Delta }𝐡^{i(4)}(\tau _i)`$ $`=`$ $`{\displaystyle \frac{\gamma \alpha _0\tau _i^2}{2!}}𝐜`$ (87)
at the most leading order. Thereafter, having combined Eqs. (75) and (83-87), all elements of $`𝐄^i(\tau _i)`$ are seen to have a prefactor of order $`\stackrel{~}{\gamma }^3`$. This makes $`𝐄^i(\tau _i)`$ small compared to $`1`$. Hence, to first order, $`𝚫L(\tau )=𝐋(\tau )𝐋_0(\tau )`$ can be found from Eq. (69) by keeping only the terms linear in $`𝐄^i(\tau _i)`$:

$`\mathrm{\Delta }𝐋(\tau )`$ $`=`$ $`e^{\frac{N1}{2}\alpha _0\tau }\left\{{\displaystyle \underset{i=1}{\overset{N}{}}}𝐇_0^i(\tau +t_0\stackrel{~}{t}_i)\right\}({\displaystyle \underset{j=1}{\overset{N}{}}}𝐄^j(\tau +t_0\stackrel{~}{t}_j){\displaystyle \underset{c=1}{\overset{s}{}}}𝐌_c𝐇_0^{i_c}(\tau _{i_c}^c)𝐇_0^{j_c}(\tau _{j_c}^c)`$ (89)
$`+{\displaystyle \underset{c=1}{\overset{s}{}}}{\displaystyle \underset{a=1}{\overset{c}{}}}𝐌_a𝐇_0^{i_a}(\tau _{i_a}^a)𝐇_0^{j_a}(\tau _{j_a}^a)\times [𝐄^{j_c}(\tau _{j_c}^c)+𝐄^{i_c}(\tau _{i_c}^c)]{\displaystyle \underset{b=c+1}{\overset{s}{}}}𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b)).`$
We now shuffle all the $`𝐄^{i_c}`$ matrices to the right, by successively exchanging the order of $`𝐄^{i_c}`$ and the next $`𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b)`$ to its right. These will commute very often, because most collisions in the sequence will involve other spheres than this $`i_c`$-th one. If they do not commute, we write

$`𝐄^{i_c}𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b)=𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b)\left[𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b)\right]^1𝐄^i𝐌_b𝐇_0^{i_b}(\tau _{i_b}^b)𝐇_0^{j_b}(\tau _{j_b}^b).`$ (90)
From that point on, we work with “$`[𝐌_b𝐇_0^{i_b}𝐇_0^{j_b}]^1`$ $`𝐄^i𝐌_b𝐇_0^{i_b}𝐇_0^{j_b}`$”, and shuffle that to the right. Repeating this, we end up with $`𝐋_0(\tau )`$ on the left side again, and find, symbolically,

$`\left[𝐋_0(\tau )\right]^1\mathrm{\Delta }𝐋(\tau )`$ $`=`$ $`{\displaystyle \underset{c=1}{\overset{s}{}}}\left\{{\displaystyle (𝐌H_0H_0)^1}\right\}\left[𝐄^{j_c}(\tau _{j_c}^c)+𝐄^{i_c}(\tau _{i_c}^c)\right]\left\{{\displaystyle 𝐌H_0H_0}\right\}`$ (92)
$`+{\displaystyle \underset{j=1}{\overset{N}{}}}\left\{{\displaystyle (𝐌H_0H_0)^1}\right\}𝐄^j(\tau +t_0\stackrel{~}{t}_j)\left\{{\displaystyle 𝐌H_0H_0}\right\}.`$
Since $`\tau `$ is of the order of a flight time of a sphere, $`s`$ is of the order of $`N`$, but a typical sphere will have suffered only a few \[a few here is roughly $`O(1)`$\] collisions. Consider a typical term in the sum over $`c`$. The $`𝐌H_0H_0`$ terms in the products are from collisions that involve the $`i_c`$-th or the $`j_c`$-th sphere, or spheres that were involved in collisions with the $`i_c`$-th or $`j_c`$-th sphere. The number of $`𝐌H_0H_0`$ will thus also be of $`O(1)`$.

The matrix $`𝐄^{i_c}`$ has non-zero elements only for the entries associated with the $`i_c`$-th sphere. The multiplication by $`𝐌H_0H_0`$ yields more non-zero entries associated with other spheres, whose number however is of $`O(1)`$. Therefore, in the whole sum of $`s=O(N)`$ terms in Eq. (92) will yield for typical entries a result of order one (and not  of order $`N`$ as the sum over $`sN`$ suggests). Since all elements of the $`𝐄`$ matrices are proportional to $`\gamma \stackrel{~}{\gamma }^2`$ as can be seen from Eqs. (75), (83) and (87), we conclude for $`tt_0=\tau =O(\tau _0)`$ that

$`𝐋(\tau )=𝐋_0(\tau )[𝐈+\stackrel{~}{\gamma }^3𝐁],`$ (93)
where the matrix $`𝐁`$ is of order 1 in $`\stackrel{~}{\gamma }`$ and order 1 in $`N`$. Note that $`𝐁`$ contains higher powers of $`\stackrel{~}{\gamma }`$ as well. Because it involves $`𝐜`$ and contributions from collisions between spheres, $`𝐁`$ is not proportional to the identity $`𝐈`$, so we cannot regard it as a simple scalar factor (in which case the exact conjugate pairing would be easy to obtain again). Equation (93) implies

$`𝚫\stackrel{~}{L}(\tau )`$ $`=`$ $`\stackrel{~}{}L(\tau )\stackrel{~}{}L_0(\tau )`$ (94)
$`=`$ $`\stackrel{~}{\gamma }^3[𝐁^\text{T}\stackrel{~}{}L_0(\tau )+\stackrel{~}{}L_0(\tau )𝐁]+\stackrel{~}{\gamma }^6𝐁^\text{T}\stackrel{~}{}L_0(\tau )𝐁.`$ (95)
We can now see, from Eqs. (93) and (95) that the differences between $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$, and between $`\stackrel{~}{}L(\tau )`$ and $`\stackrel{~}{}L_0(\tau )`$ are small, as anticipated on page IV C, by a relative order $`\stackrel{~}{\gamma }^3`$. Therefore the logarithm of the eigenvalues of $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$ also differ by a term of order $`\stackrel{~}{\gamma }^3`$ in an absolute sense. If we now divide the logarithms of these eigenvalues by the time $`\tau `$, we see that the finite time (for time $`\tau `$) Lyapunov exponents, calculated from $`\stackrel{~}{}L_0`$ and from $`\stackrel{~}{}L`$ (which we denote as $`\lambda _i^{(0)}(\tau )`$ and $`\lambda _i(\tau )`$ respectively, for $`i=1,2\mathrm{}6N`$) differ by a term $`O(\stackrel{~}{\gamma }^3/\tau )=O(\gamma \stackrel{~}{\gamma }^2)`$. Using the fact that that $`\lambda _i^{(0)}(\tau )`$s satisfy the conjugate pairing rule, $`\lambda _i(\tau )`$s too, will satisfy it up to $`O(\gamma \stackrel{~}{\gamma }^2)`$. From the condition that

$`{\displaystyle \underset{i=1}{\overset{6N}{}}}\lambda _i^{(0)}(\tau )={\displaystyle \underset{i=1}{\overset{6N}{}}}\lambda _i(\tau )=3N\alpha _0,`$ (96)
\[which can be easily verified from Eqs. (42) and (58)\], we see that on average over the pairs, the deviation from CPR is zero.

We make one further observation at this stage. The Lyapunov exponents (even the finite time once) are invariant under $`\gamma \gamma `$, so in a power series expansion of the individual Lyapunov exponents in $`\stackrel{~}{\gamma }`$ , the odd powers vanish. We therefore conclude that the logarithm of the eigenvalues of $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$ must differ by a term of order $`\stackrel{~}{\gamma }^4`$, and hence conjugate pairing of $`\lambda _i(\tau )`$s will be valid up to the correction of the form $`\gamma \stackrel{~}{\gamma }^3`$.

To extend these ideas for large $`(tt_0)`$ \[and finally for $`t\mathrm{})`$\], we proceed in the following way. Notice that the matrices $`\stackrel{~}{}L(tt_0)`$ and $`\stackrel{~}{}L_0(tt_0)`$ are positive definite and symmetric. This allows us to express the two matrices in the form of $`\stackrel{~}{}L_0(tt_0)=\mathrm{exp}(𝐀_0)`$ and $`\stackrel{~}{}L(tt_0)=\mathrm{exp}(𝐀)`$, where for large $`(tt_0)`$, both the eigenvalues of $`𝐀_0`$ and $`𝐀`$ must behave $`(tt_0)`$. In these terms, the difference between the Lyapunov exponents for $`\stackrel{~}{}L(tt_0)`$ and $`\stackrel{~}{}L_0(tt_0)`$ is related to $`𝐀𝐀_0`$. Since the difference between $`\stackrel{~}{}L(tt_0)`$ and $`\stackrel{~}{}L_0(tt_0)`$ has an explicit prefactor of $`\stackrel{~}{\gamma }^3`$, so will $`𝐀𝐀_0`$. Using the symmetry argument that the Lyapunov exponents have to be even functions of $`\gamma `$, we obtain

$`\lambda _i+\lambda _{6Ni+1}=\alpha _0+O(\gamma \stackrel{~}{\gamma }^3).i=1,\mathrm{}6N`$ (97)
To explicitly extend the formalism developed in Eqs. (69-95) to large $`tt_0`$ and thereby obtain a relation between $`\lambda _i`$s and $`\lambda _i^{(0)}`$s, we need to concatenate a lot of $`𝐋(\tau )`$s. These do not commute with each other, nor do they commute with $`𝐁`$s in general. This prevents an explicit demonstration of how the deviation is built up. For the largest and the most negative Lyapunov exponent, it has been possible to show that they pair to $`\alpha _0`$ plus corrections of $`O(\gamma \stackrel{~}{\gamma }^3)`$ by means of a kinetic theory approach , based on the independence of subsequent collisions of a sphere. One expects that in subsequent time-intervals of $`O(\tau _0)`$, the $`𝐋(\tau )`$ matrices are not qualitatively much different from each other. Therefore, we expect that the coefficient of the $`O(\gamma \stackrel{~}{\gamma }^3)`$ term in Eq. (97), to be of the same order as that for a flight time $`\tau =O(\tau _0)`$, i.e. of the order of $`𝐁`$, which is $`O(1)`$.

We now return to the discussion of an approximate CPR for the SLLOD equations of motion with an iso-kinetic Gaussian thermostat in between collisions. As discussed in Sec. III, in the non-equilibrium steady state and in the thermodynamic limit, the coefficient of dynamical friction representing the iso-kinetic Gaussian thermostat for the SLLOD dynamics fluctuates with $`1/\sqrt{N}`$ fluctuations around the fixed value $`\alpha _0`$. We would therefore expect that the Lyapunov spectrum for the SLLOD dynamics with an iso-kinetic Gaussian thermostat are given by that of $`𝐋(tt_0)`$ plus  terms of $`O(1/\sqrt{N})`$. Consequently, the approximate conjugate pairing of Eq. (97) can be extended to this thermostat once we neglect the $`O(1/\sqrt{N})`$ terms in the sum of the corresponding $`\lambda _i`$ and $`\lambda _{6Ni+1}`$.

## V Discussion

In this paper, we started with a collection of $`N`$ hard spheres, each with unit mass and arbitrary radius. Next, we argued how the coefficient of friction representing the iso-kinetic Gaussian thermostat in between collisions, fluctuates around a fixed value $`\alpha _0`$ with $`1/\sqrt{N}`$ fluctuations in the NESS. Using the properties of the transformation matrices for the infinitesimal phase space element $`\delta 𝚪`$, we then showed that the CPR is expected to be violated, at the most, at $`O(\gamma ^4)`$ for constant coefficient of friction $`\alpha _0`$. The source of the CPR violation is basically the $`\alpha _0\gamma \mathrm{𝐂𝐑}`$ term in Eq. (7), as it is clearly seen that CPR would have been exact in the absence of this term. Finally, we extended that result to the case when the coefficient of friction represents an iso-kinetic Gaussian thermostat. Note that to obtain the deviations from CPR for the constant $`\alpha _0`$ thermostat, the number of particles did not matter, so for that case the result holds for any number of particles, but with an iso-kinetic Gaussian thermostat between collisions, we need to take $`N`$ to infinity to make the connection to the constant $`\alpha _0`$ thermostat.

In addition, as mentioned in the introduction of this paper, the condition that the mass of each sphere be unity, is not necessary. If we assume that the mass of the $`i`$-th sphere is $`m_i`$, then one can obtain the same equations of motion, Eq. (1-2) and (7-12) in terms of the primed quantities defined by $`𝐫_i^{}=\sqrt{m_i}𝐫`$, $`𝐯_i^{}=\sqrt{m_i}𝐯`$, $`𝐩_i^{}=𝐩_i/\sqrt{m_i}`$, $`𝐑^{}=(𝐫_1^{},𝐫_2^{},\mathrm{},𝐫_N^{})`$, $`𝐕^{}=(𝐯_1^{},𝐯_2^{},\mathrm{},𝐯_N^{})`$ and $`\widehat{𝐍}_{ij}^{}`$, defined as

$`\widehat{𝐍}_{ij}^{}_{}{}^{}l=\left(\delta _{l,i}\sqrt{{\displaystyle \frac{m_j}{m_i+m_j}}}\delta _{l,j}\sqrt{{\displaystyle \frac{m_i}{m_i+m_j}}}\right)\widehat{𝐧}_{ij}`$

$`(l=1,2,..,N)`$. It is then straightforward, albeit lengthy and laborious, to see that our entire analysis goes through in terms of these primed variables, once we use $`𝐔^{}`$ in Eq. (A22) with

$`𝐔_{ij}^{}`$ $`=`$ $`𝐔_{ji}^{}=[m_i+m_j]^{1/2}𝐈,`$ (98)
$`𝐔_{ii}^{}`$ $`=`$ $`{\displaystyle \frac{m_j}{m_i}}𝐔_{jj}^{}=\sqrt{{\displaystyle \frac{m_j}{m_i}}}𝐔_{ij}^{}.`$ (99)
However, we must note that even though the analysis presented in Secs. IV is not affected when the masses and the radii of the spheres are not necessarily the same, one should not allow extreme variations in masses and the radii of the spheres. For large variation of the masses, the system may phase separate into phases in which the mean flight times are different, which will invalidate the use of the Boltzmann equation in Ref.  (the Boltzmann equation has been used to show the approach of $`\alpha `$ towards the constant value $`\alpha _0`$ in Ref. ). Furthermore, the use of a typical mean free time in Sec. IV.B will not be possible for large variation in the radii of the spheres. Our analysis in this paper therefore, holds for somewhat limited variation of the masses and the radii of the spheres.

Another possible generalization of our analysis can be carried out for the case when the gas particles interact with each other by means of a short-ranged, repulsive, inter-particle potential (an attractive potential may cause bound states) for a constant multiplier thermostat. The dynamics of the particles can be again decomposed into flights and “collisions” at low density. While the transformation of the infinitesimal phase space volume element over flights will again be determined by the $`𝐇`$ matrices as in Eqs. (37-42) (and thereby have the same properties as above), the corresponding $`𝐌_{ij}`$ matrices will not necessarily have a similar form and properties as presented in Eqs. (58-59). For a given average value of the peculiar kinetic energy, the constant multiplier thermostat has to be chosen carefully (see the last paragraph of Sec. III). Having denoted the constant coefficient of friction again by $`\alpha _0`$, where $`\alpha _0\gamma ^2`$, is easy to see that the resulting matrix $`𝐌_{ij}`$ can be decomposed into a sum of two matrices, $`𝐌_{ij}^{(0)}`$ and $`𝚫M_{ij}`$, where $`𝐌_{ij}^{(0)}`$ is exactly $`\mu `$-symplectic but $`𝚫M_{ij}`$ is not. If the term $`\alpha _0\gamma \mathrm{𝐂𝐑}`$ in Eq. (7) (which is the culprit for the violations of CPR in the first place) had been absent, then $`𝐌_{ij}`$ would have been $`\mu `$-symplectic  and $`𝚫M_{ij}`$ would have been zero, but with Eq. (7) in its present form, it is easily seen that $`𝚫M_{ij}`$ is $`\alpha \gamma \stackrel{~}{\gamma }^3`$ and the proportionality constant depends on the ratio of the time the two colliding particles spend in contact with each other to the time of a flight. This ratio is very small at low density of the gas and therefore the elements of $`𝚫M_{ij}`$ are small in comparison with $`𝐌_{ij}^{(0)}`$. Since the $`𝐇`$ matrices have a similar property, one can repeat the analysis of this paper for such inter-particle short-range repulsive potentials using $`𝐇_0`$ and $`𝐌_{ij}^{(0)}`$ as the reference matrices and arrive at the same conclusion as Eq. (97).

Moreover, on the basis of the fact that at low densities and short-ranged repulsive interparticle potentials (that do not give rise to bound states), the constant multiplier thermostat and an iso-kinetic Gaussian thermostat are equivalent in the thermodynamic limit , one would expect that under such conditions, the Lyapunov exponent spectrum for the Gaussian thermostat would be practically the same as that of the corresponding constant multiplier thermostat, where the constant multiplier is once again proportional to $`\gamma ^2`$. Therefore, in view of the discussion in the previous paragraph, we can also expect $`O(\stackrel{~}{\gamma }^4)`$ deviations from the CPR for a Gaussian thermostatted SLLOD equations of motion in the thermodynamic limit, where the particles interact with each other by means of a short-ranged repulsive potential.

Finally, we note here that verifications of our theory by means of computer simulation remains a challenge. The work for simulations is in progress at present, although our preliminary experience suggests that to retrieve the $`\gamma \stackrel{~}{\gamma }^3`$ scaling of the deviation of the sum of the pairs of Lyapunov exponents from $`\alpha _0`$ is not an easy task.

## Acknowledgment

The authors wish to thank Prof. J. R. Dorfman and Prof. H. van Beijeren for many useful and motivating discussions regarding this subject. D. P. was supported by the research grants from “Fundamenteel Onderzoek der Materie (FOM)” and that of Prof. J. R. Dorfman, NSF-PHY-9600428, and wishes to acknowledge the kind hospitality of Prof. H. van Beijeren at University of Utrecht during this and related work. R. v. Z. was supported by the research grant of Prof. H. van Beijeren, by FOM, SMC and NWO Priority Program Non-Linear Systems, which are financially supported by the ”Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)”, and by a grant from the Natural Sciences and Engineering Research Council of Canada.

##

In this appendix, we will derive the transformation $`𝐌_{ij}`$ of tangent vectors at a collision. The essential element is that the collision does not happen at the same time on reference and on the adjacent trajectory. To understand the origin of this time lag between the binary collisions between the $`i`$-th and the $`j`$-th sphere on the reference and the adjacent trajectories, in Fig. 1, we have depicted an exaggerated schematic picture of the collisions taking place in the $`6N`$-dimensional phase space on the reference and the adjacent trajectories at points A and C respectively. The points B and D show the corresponding positions of the adjacent and the reference points respectively when the binary collisions at A and C are taking place. Thus, the pre-collisional separation between the reference and the adjacent points is $`\stackrel{}{\text{AB}}=\delta 𝚪_{}`$, while the post-collisional separation is $`\stackrel{}{\text{DC}}=\delta 𝚪_+`$. Using that $`|𝐫_j𝐫_i|`$ and $`|𝐫_j+\delta 𝐫_j+𝐯_j\delta \tau 𝐫_i\delta 𝐫_i𝐯_i\delta \tau |`$ both have to equal $`a_i+a_j`$, the time lag $`\delta \tau `$ between the two collisions at A and C, therefore, can be easily expressed as

$`\delta \tau ={\displaystyle \frac{(\delta 𝐫_j\delta 𝐫_i)\widehat{𝐧}_{ij}}{(𝐯_j𝐯_i)\widehat{𝐧}_{ij}}}={\displaystyle \frac{\delta 𝐑_{}\widehat{𝐍}_{ij}}{𝐕_{}\widehat{𝐍}_{ij}}}.`$ (1)
To obtain the expression of $`\delta 𝚪_+`$, we first express the transformation of $`(𝐑,𝐕)`$ in Eq. (12) over a binary collision between the $`i`$-th and the $`j`$-th sphere in a matrix form

$`𝚪_+=𝐐(𝚪_{})=\left[\begin{array}{c}𝐑_{}\\ 𝐕_{}2(𝐕_{}\widehat{𝐍}_{ij})\widehat{𝐍}_{ij}\end{array}\right],`$ (4)
Next we notice that

$`\delta 𝚪_+={\displaystyle \frac{𝐐}{𝚪_{}}}\delta 𝚪_{}^{}\dot{𝚪}_+\delta \tau ,`$ (5)
where, $`\delta 𝚪_{}^{}`$ is the infinitesimal phase space separation between the two trajectories at A and C is

$`\delta 𝚪_{}^{}=\delta 𝚪_{}+\dot{𝚪}_{}\delta \tau `$ (6)
and $`\dot{𝚪}_\pm `$ describes the equations of motion, Eq. (7), right after(before) the collision at A, i.e.,

$`\dot{𝚪}_\pm =\left[\begin{array}{c}\dot{𝐑}_\pm \\ \dot{𝐕}_\pm \end{array}\right]=\left[\begin{array}{c}𝐕_\pm \\ \alpha _0\gamma \mathrm{𝐂𝐑}_\pm \alpha _0𝐕_\pm \end{array}\right].`$ (11)
Having calculated the quantity $`{\displaystyle \frac{𝐐}{𝚪_{}}}`$ from Eq. (4),

$`{\displaystyle \frac{𝐐}{𝚪_{}}}=\left[\begin{array}{cc}𝐈& \mathrm{𝟎}\\ \mathrm{\hspace{0.17em}2}𝐕_{}\left[{\displaystyle \frac{\widehat{𝐍}_{ij}}{𝐑_{}}}\widehat{𝐍}_{ij}+\widehat{𝐍}_{ij}{\displaystyle \frac{\widehat{𝐍}_{ij}}{𝐑_{}}}\right]& 𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij}\end{array}\right],`$ (A14)
where each entry of the matrix on the r.h.s. of Eq. (A14) is a $`3N\times 3N`$ matrix, the expressions for $`\delta 𝐑_+`$ and $`\delta 𝐕_+`$ can be obtained from Eq. (5) as

$`\delta 𝐑_+=\delta 𝐑_{}\mathrm{\hspace{0.17em}2}(\delta 𝐑_{}\widehat{𝐍}_{ij})\widehat{𝐍}_{ij},`$ (A15)
and

$`\delta 𝐕_+=𝐀\delta 𝐑_{}^{}+(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})\delta 𝐕_{}+(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})[\alpha _0\gamma \mathrm{𝐂𝐑}_{}\alpha _0𝐕_{}]\delta \tau [\alpha _0\gamma \mathrm{𝐂𝐑}_+\alpha _0𝐕_+]\delta \tau ,`$ (A16)
where $`𝐀=2𝐕_{}\left[{\displaystyle \frac{\widehat{𝐍}_{ij}}{𝐑_{}}}\widehat{𝐍}_{ij}+\widehat{𝐍}_{ij}{\displaystyle \frac{\widehat{𝐍}_{ij}}{𝐑_{}}}\right]`$ and $`\delta 𝐑_{}^{}=(\delta 𝐑_{}+𝐕_{}\delta \tau )`$. At this point, we use Eq. (12) and obtain

$`(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})[\alpha _0\gamma \mathrm{𝐂𝐑}_{}\alpha _0𝐕_{}]\delta \tau [\alpha _0\gamma \mathrm{𝐂𝐑}_+\alpha _0𝐕_+]\delta \tau =\mathrm{\hspace{0.17em}2}\alpha _0\gamma (\widehat{𝐍}_{ij}\mathrm{𝐂𝐑}_{})\widehat{𝐍}_{ij}\delta \tau .`$ (A17)
Following Appendix B of , the term $`𝐀\delta 𝐑_{}^{}`$ can be expressed as

$$𝐀\delta 𝐑_{}^{}=2[(𝐕_{}\delta 𝐍_{ij})\widehat{𝐍}_{ij}+(𝐕_{}\widehat{𝐍}_{ij})\delta 𝐍_{ij}],$$
(A18)
where

$`\delta 𝐍_{ij}={\displaystyle \frac{1}{\sqrt{2}}}(\mathrm{𝟎},\mathbf{\hspace{0.17em}0},..,\delta 𝐧_{ij},\mathrm{}\delta 𝐧_{ij},\mathrm{}\mathrm{𝟎}),`$ (A19)
satisfying $`\widehat{𝐧}_{ij}\delta 𝐧_{ij}=0`$. This orthogonality condition between $`\widehat{𝐧}_{ij}`$ and $`\delta 𝐧_{ij}`$ also implies that $`\widehat{𝐍}_{ij}\delta 𝐍_{ij}=\mathrm{\hspace{0.17em}0}`$. To obtain an expression for $`\delta 𝐧_{ij}`$, we need to take a look at Figs. 2 and 3. Figure 2 describes, in the laboratory frame, the binary collision process between the $`i`$-th and the $`j`$-th sphere on the reference and adjacent trajectories; the thick-lined spheres are on the reference trajectory whereas the thin-lined spheres are on the adjacent trajectory. Figure 3 describes the same binary collision process in the reference frame of the $`i`$-th sphere (with center C). In Fig. 2, the thick-lined $`j`$-th sphere (with center D) on the left depicts the collision situation on the reference trajectory and the thin-lined $`j`$-th sphere (with center E) on the left depicts the collision situation on the adjacent trajectory.

Clearly, in Fig. 3, the infinitesimal vector $`\stackrel{}{\text{DE}}`$ is given by

$`\delta 𝐫_{ij}^{}=\delta 𝐫_j\delta 𝐫_i+(𝐯_j𝐯_i)\delta \tau `$ (A20)
and since the lengths of both the lines CD and CE are $`a_i+a_j`$ ($`a_i`$ and $`a_j`$ are the radii of the $`i`$-th and the $`j`$-th sphere respectively), we have

$`\delta 𝐧_{ij}={\displaystyle \frac{1}{a_i+a_j}}\delta 𝐫_{ij}^{}.`$ (A21)
Let us define a $`3N\times 3N`$ matrix $`𝐔`$ composed of $`N\times N`$ blocks of $`3\times 3`$ matrices, such that, in terms of the block indices the only non-zero entries of $`𝐔`$ are

$`𝐔_{ii}=𝐔_{ij}=𝐔_{ji}=𝐔_{jj}=𝐈,`$ (A22)
where $`𝐈`$ is the $`3\times 3`$ unit matrix. One can now write, using Eqs. (A19-A22) and Eq. (1) that

$`\delta 𝐍_{ij}={\displaystyle \frac{1}{\sqrt{2}(a_i+a_j)}}[𝐔\delta 𝐑_{}+𝐔\dot{𝚪}_{}\delta \tau ]={\displaystyle \frac{1}{\sqrt{2}(a_i+a_j)}}\left[𝐔{\displaystyle \frac{(𝐔𝐕_{})\widehat{𝐍}_{ij}}{𝐕_{}\widehat{𝐍}_{ij}}}\right]\delta 𝐑_{}.`$ (A23)
From Eqs. (A18) and Eq. (A23), we finally have

$`𝐀\delta 𝐑_{}^{}`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{a_i+a_j}}\left[(𝐔𝐕_{})\widehat{𝐍}_{ij}+\widehat{𝐍}_{ij}(𝐔𝐕_{})(\widehat{𝐍}_{ij}𝐕_{})𝐔+{\displaystyle \frac{𝐕_{}𝐔𝐕_{}}{𝐕_{}\widehat{𝐍}_{ij}}}\widehat{𝐍}_{ij}\widehat{𝐍}_{ij}\right]\delta 𝐑_{}`$ (A24)
$`=`$ $`(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})𝐖\delta 𝐑_{}.`$ (A25)
where $`𝐖`$ is a $`3N\times 3N`$ symmetric matrix. Finally, using Eqs. (1), (A15-A17) and (A25), the expression for $`𝐌_{ij}`$ can be obtained as

$`𝐌_{ij}=(𝐈2\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})\left[\begin{array}{cc}𝐈& \mathrm{𝟎}\\ 𝐑& 𝐈\end{array}\right],`$ (A28)
where

$$𝐑=𝐖2\alpha _0\gamma \frac{\widehat{𝐍}_{ij}\mathrm{𝐂𝐑}_{}}{𝐕_{}\widehat{𝐍}_{ij}}\widehat{𝐍}_{ij}\widehat{𝐍}_{ij}.$$
(A29)


text:
# A Fabry-Perot Study of the Scd I galaxy NGC 5457

## 1. Introduction

The study of the kinematics of spiral galaxies is an active area of astronomical research. The construction of Fabry-Perot interferometers have greatly increased the kinematical knowledge of these objects. In this contribution, we present preliminary results from a Fabry-Perot study of NGC 5457 (M101).

## 2. Data and Reduction

The observations were carried out during the night of July 23-24, 2001, with the UNAM Scanning Fabry-Perot interferometer (SFPI) PUMA attached to the f/7.9 Ritchey-Chretien focus of the 2.1m telescope at the Observatorio Astronómico Nacional at San Pedro Mártir, B.C., México. The main characteristics of PUMA can be found in Rosado et al. (1995). The data reduction was performed using the ADHOC package (Boulesteix 1993). The reduction procedure followed a “standard” scheme and the details will be published elsewhere (Puerari et al., in preparation).

## 3. Results

One of the main results obtained from Fabry-Perot data is the radial velocity field. By using this field, and the geometrical parameters of the galaxy (inclination and position angles, kinematical center) as well as the recession velocity, the rotation curve can be calculated. For NGC 5457, we have derived the rotation curve presented in Fig. 1 (left panels). This curve is in agreement with that one of Comte et al. (1979, see their Fig. 6), but our results present less dispersion, due to the better resolution and signal/noise ratio of our observations. Assuming a distance of 7.2 Mpc (Sandage & Tammann 1974), we have calculated a total mass of 9.8$`\times `$ 1010 M inside a radius of 4$`.^{}`$8 (or 10 kpc), also in agreement with Comte et al. (1979).

We have also obtained the velocity dispersion of the HII region population. We fitted a gaussian to each HII region velocity profile and corrected the final value by the instrumental, thermal and intrinsic broadenings. The distribution of velocity dispersion values presents a normal behavior (see Fig. 1, right panel). The mean value of this distribution is 30 km sec-1 (i.e., supersonic dispersion); this result must be checked with Fabry-Perot studies at other lines (eg., \[SII\]).

## References

Boulesteix, J. 1993, “ADHOC Reference Manual”, Publications de l’Observatoire de Marseille (http://alpha2.cnrs-mrs.fr/adhoc/)

Comte, G., Monnet, G., Rosado, M. 1979, A&A, 72, 73

Rosado, M. et al. 1995, Rev. Mex. Astron. y Astrof., 3, 263

Sandage, A., Tammann, G.A. 1974, ApJ, 194, 223


text:
# Ξ⁻ and Ξ̄+ Production in Central Pb+Pb Collisions at 158 GeV/𝑐 per Nucleon

## 1 Introduction

High energy heavy–ion physics provides a unique opportunity to study nuclear matter under extreme conditions of temperature and pressure. Quantum Chromodynamics (QCD) predicts that at high energy density, nuclear matter will melt into a deconfined state of partonic matter known as the Quark–Gluon Plasma (QGP). Experimental efforts to detect this novel phase rely on measurements of final state particles emitted from the fireball created by colliding nuclear matter. There are three key stages in the collision process in which a QGP is formed: the initial partonic interactions which occur as two nuclei interpenetrate, the subsequent partonic interactions in a thermalised QGP and finally hadronic interactions once the temperature of the system has dropped below the critical temperature. In order to probe the QGP, signals must be sought which are produced dominantly in the equilibrated deconfined state and survive the hadronisation process. (For a general overview, we refer to the proceedings of the Quark Matter ’99 conference .)

For almost 20 years, the production of strangeness has been proposed as a sensitive signal to QGP formation . (For a recent overview see the proceedings of the Strangeness 2000 conference .) Strange particles are of particular interest in hadronic collisions since they carry a new quantum number not present in the colliding nucleons or nuclei. In a QGP, $`s\overline{s}`$ pairs are readily created because of the newly available gluonic degrees of freedom which contribute approximately 80$`\%`$ to the total strangeness produced . Strangeness in a partonic system where chiral symmetry is partially restored, is also expected to be readily produced due to the lower energy threshold for $`s\overline{s}`$ quark pair production. At finite baryon density, where the chemical potential associated with the production of further light quarks is raised, the production of strangeness may be further increased. These considerations lead to the expectation that the QGP phase should be characterised by a high strange (anti)quark density . By virtue of their strange quark content, multi–strange baryons and anti–baryons are expected to be a more sensitive probe to this phase than hadrons which contain only one strange valence quark. The production of multi–strange baryons, and particularly multi–strange anti–baryons is otherwise suppressed by high hadronic energy thresholds as well as by long timescales for multi–step processes. The NA49  and WA97  collaborations have previously reported an enhancement of multi–strange hyperons from central Pb+Pb collisions at 158 GeV/$`c`$ per nucleon at midrapidity compared to a convolution of nucleon-nucleon collisions. Total yields of strange hadrons provide additional information, since they are an essential ingredient in thermal models . Remarkably, these models have described a wide range of experimental data with much success . Notable exceptions are in the multi–strange sector . This paper reports on the rapidity distributions of $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ measured by the NA49 collaboration over a wider range than has previously been published.

## 2 The NA49 Experiment and Global Analysis

The NA49 large acceptance hadron spectrometer has been used to record the charged final state particles produced in central Pb+Pb collisions at 158 GeV per nucleon at the CERN SPS. The main tracking devices are four large volume time projection chambers (TPC) which are capable of detecting 80$`\%`$ of some 1500 charged particles created in a central Pb+Pb collision. Two TPCs, VT1 and VT2, are located in strong magnetic fields of 1.5 and 1.1 T respectively, which allows charge and momentum to be assigned to the detected particles. A further two TPCs are situated downstream from the magnets in a field free region. The results presented here were analysed with a global tracking scheme , which combines track segments that belong to the same physical particle but are detected in different TPCs. Generally, global tracks have better momentum resolution. Utilising a global tracking scheme increases detector acceptance and improves the track reconstruction efficiency in the TPCs, because tracks measured in the lower density regions far from the interaction vertex can be used to identify track segments in high density regions. Full details of the NA49 detector set–up and performance of tracking are available in .

The results contained in this paper were obtained from the analysis of 400,000 central Pb+Pb events. The events were selected according to the energy deposited in a forward calorimeter and correspond to the 10$`\%`$ most central collisions with an average number of participants $`N_{part}`$ $``$ 335.

## 3 $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ from Central Pb+Pb Collisions

The doubly–strange $`\mathrm{\Xi }^{}`$ hyperon decays via the channel $`\mathrm{\Xi }^{}\mathrm{\Lambda }+\pi ^{}`$ ($`100\%`$ branching ratio) with the subsequent decay $`\mathrm{\Lambda }p+\pi ^{}`$ ($`63.9\%`$ branching ratio). Neutral particle decay candidates were formed by considering all possible pairs of positively and negatively charged tracks reconstructed in the TPCs. Each pair of tracks was extrapolated from their most upstream point back towards the target. If the distance of closest approach between the two tracks was less than 1 cm at any point along their trajectory the pair was saved as a candidate $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$). A large fraction of the combinatorial background arising from random pairs of tracks was eliminated by requiring that the tracks were separated by at least 2 cm at the target plane and that the decay vertex was at least 100 cm downstream of the target. If the invariant mass of the candidates fell within $`\pm 15`$ MeV/$`c^2`$ of the $`\mathrm{\Lambda }`$ ($`\overline{\mathrm{\Lambda }}`$) mass, they were selected as candidate daughters of $`\mathrm{\Xi }^{}`$($`\overline{\mathrm{\Xi }}`$+) decays. Each $`\mathrm{\Lambda }`$ candidate was paired with each negatively charged track in the event, and both were extrapolated back toward the target. As in the case of reconstructing the neutral decay vertices, if the distance of closest approach between trajectories of the neutral $`\mathrm{\Lambda }`$ candidate and the negatively charged track fell within 1 cm, the pair was saved as a possible $`\mathrm{\Xi }^{}`$ candidate. Equivalently, each $`\overline{\mathrm{\Lambda }}`$ candidate was paired with each positive track in the event. Once again, straight-forward geometrical cuts were applied to eliminate most of the combinatorial background. In this case, the most important cuts were that the trajectories of the charged decay daughters were separated by between 2.5 cm and 7 cm from the primary vertex in the target plane depending on the position of the decay vertex. An additional requirement was that the decay vertex was at least 35 cm downstream of the target and that the $`\mathrm{\Xi }^{}`$($`\overline{\mathrm{\Xi }}`$+) candidate pointed back to the primary vertex. The resulting invariant mass distributions of $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ candidates are shown in figure 1.

The signal was obtained by selecting those particles with invariant mass within $`\pm 20`$ MeV/$`c^2`$ of the $`\mathrm{\Xi }^{}`$ mass and then subtracting background counts from regions outside the peak. Specifically, the background regions were defined between the ranges 1267 $`<M<`$ 1284 MeV/$`c^2`$ and 1358 $`<M<`$ 1381 MeV/$`c^2`$. The background subtracted data set used consists of approximately 4800 $`\mathrm{\Xi }^{}`$ and 900 $`\overline{\mathrm{\Xi }}`$+ reconstructed particles .

Those $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ which were excluded by the geometric selection criteria, or not found due to reconstruction losses, were later accounted for by applying a correction for the overall reconstruction inefficiency.

### 3.1 Transverse Mass Distributions

The data were evaluated in bins of rapidity, transverse momentum and lifetime. Corrections for geometrical acceptance, branching ratio and reconstruction efficiency were then applied on a bin–by–bin basis. The GEANT simulation package has been adapted for NA49 and was used for calculating the geometrical acceptance. For the determination of the reconstruction efficiency, a single Monte Carlo $`\mathrm{\Xi }^{}`$ (or $`\overline{\mathrm{\Xi }}`$+) generated within the spectrometer acceptance was embedded into the raw data structures in which no real $`\mathrm{\Xi }^{}`$ (or $`\overline{\mathrm{\Xi }}`$+) was detected. These events were processed with the reconstruction software which tried to reconstruct the embedded particle. The reconstruction efficiency of the embedded $`\mathrm{\Xi }^{}`$ (or $`\overline{\mathrm{\Xi }}`$+) was determined statistically for each bin of rapidity, transverse momentum and lifetime. The kinematic phase space populated by detectable $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ hyperons was found to be $`0.6<p_T<2.8`$ GeV/$`c`$ in transverse momentum and $`1.7<y<4.5`$ in rapidity. Note that the rapidity range covers both sides of midrapidity ($`y_{\mathrm{cm}}=2.9`$). Fully corrected transverse mass ($`m_T=\sqrt{p_{T}^{}{}_{}{}^{2}+m_{0}^{}{}_{}{}^{2}}`$) distributions are shown in figure 2 and are fitted with the function given in equation (1) to obtain the inverse slope parameters, $`T`$.

$$\frac{1}{m_T}\frac{dN}{dm_T}\mathrm{exp}\left(\frac{m_T}{T}\right)$$
(1)
The measured transverse mass spectra span the detectable rapidity range in each $`p_T`$ bin. The fitted temperature parameters are found to be 267 $`\pm `$ 9 MeV and 296 $`\pm `$ 17 MeV for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ hyperons respectively.

### 3.2 Rapidity Distributions, 4$`\pi `$ Yields and Particle Ratios

The large acceptance coverage of the NA49 spectrometer allows the experimental determination of fully integrated particle yields. However, physical measurements are not possible for all values of transverse momentum because of finite geometrical acceptance. Consequently, the rapidity spectra for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ are formed by extrapolating the transverse momentum distributions over the full range by means of eq. (1) using the fitted slope parameters, and assuming factorisation in $`p_T`$ and $`y`$. The resulting distributions are shown in figure 3. The extrapolation factors are 1.51 and 1.45 for the $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ distributions, respectively. At midrapidity, the yields per event per unit of rapidity (or central rapidity densities, $`dN/dy`$) are found to be $`1.49\pm 0.08`$ and $`0.33\pm 0.04`$ for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ hyperons respectively.

Integration of the measured part of the rapidity distribution and taking into account the reflected data, averaging where necessary, provides a lower limit for the cascade multiplicities of $`3.76\pm 0.12`$ and $`0.72\pm 0.05`$ particles per event for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ respectively. In order to estimate the 4$`\pi `$ yield, a fit to a double gaussian has been made in the case of the $`\mathrm{\Xi }^{}`$ rapidity spectrum. Each gaussian has a width $`\sigma =0.75\pm 0.12`$ centred at $`y_{cm}\pm 0.66`$. The integral of this fit over the full rapidity range provides an estimated 4$`\pi `$ yield of $`4.12\pm 0.20`$ $`\mathrm{\Xi }^{}`$ per event. A single Gaussian centred at midrapidity with a width of $`\sigma =0.87\pm 0.07`$ is shown in figure 3 for the $`\overline{\mathrm{\Xi }}`$+ distribution. Integrating this fit over the full rapidity range gives an estimated 4$`\pi `$ yield of $`0.77\pm 0.04`$ $`\overline{\mathrm{\Xi }}`$+ per event. The $`\overline{\mathrm{\Xi }}`$$`{}_{}{}^{+}/\mathrm{\Xi }^{}`$ ratio, $`R`$, can be calculated over the full measured rapidity range. At midrapidity, $`R`$($`y`$=$`y_{\mathrm{cm}}`$) is found to be $`0.22\pm 0.03`$ in good agreement with our previous publication  and with WA97 . Note that in reference , the event selection was based on the 5$`\%`$ most central collisions (10$`\%`$ in this paper), where $`N_{part}`$ was approximately 362 compared to 335 in this analysis. Additionally, the large acceptance of the NA49 spectrometer allows a determination of the 4$`\pi `$ integrated ratio, which is found to be $`R`$(4$`\pi `$) $`=0.17\pm 0.01`$.

### 3.3 Lifetime Distributions

In addition to dividing the data into $`y`$ and $`p_T`$ bins, the data were also evaluated in discrete bins of lifetime, $`\tau `$, allowing lifetime distributions to be measured. Corrections for geometrical acceptance and reconstruction efficiency were applied in each of the lifetime bins and the resulting corrected distributions are shown in figure 4. The distributions were fitted with exponentials to determine the mean lifetimes, which were found to be $`5.2\pm 0.2`$ cm/$`c`$ and $`5.0\pm 0.3`$ cm/$`c`$ for $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ hyperons respectively. These values are consistent with the data book mean lifetime value of $`\tau _0=`$ 4.91 cm/$`c`$ .

### 3.4 Discussion and Conclusions

This new analysis uses a completely revised reconstruction procedure and covers a wide range in rapidity. The resulting total integrated yields of $`\mathrm{\Xi }^{}`$ are 30$`\%`$ lower than previous NA49 estimates . This is due in part to the different centrality selection. The current measurement benefits from a smaller error in extrapolation since we have been able to measure for the first time, the rapidity distributions of both $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ over a wider range than previously. The central rapidity densities are in agreement with those found by WA97 . Errors on the new results are statistical only, with the systematic errors estimated to be about 15$`\%`$ derived from differences between simulation and data.

In summary, $`\mathrm{\Xi }^{}`$ and $`\overline{\mathrm{\Xi }}`$+ hyperons have been measured at 158 GeV/$`c`$ per nucleon using the NA49 large acceptance spectrometer. For the $`\mathrm{\Xi }^{}`$ hyperons, the slope parameter determined from a fit to the transverse mass spectrum is 267 $`\pm `$ 9 MeV, the central rapidity density is found to be $`1.49\pm 0.08`$ and the total yield is found to be $`4.51\pm 0.20`$ particles per event integrated over all phase space. The $`\overline{\mathrm{\Xi }}`$+ hyperons are found to have a slope parameter of 296 $`\pm `$ 17 MeV, with a central rapidity density of $`0.33\pm 0.04`$ and a fully integrated yield of $`0.77\pm 0.04`$ particles per event. At midrapidity, the $`\overline{\mathrm{\Xi }}`$$`{}_{}{}^{+}/\mathrm{\Xi }^{}`$ ratio is found to be $`R`$($`y`$=$`y_{\mathrm{cm}}`$) $`=0.22\pm 0.03`$. A full phase space integrated $`\overline{\mathrm{\Xi }}`$$`{}_{}{}^{+}/\mathrm{\Xi }^{}`$ ratio of $`R`$(4$`\pi `$) $`=0.19\pm 0.01`$ is obtained.

### 3.5 Acknowledgements

This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the US Department of Energy (DE-ACO3-76SFOOO98 and DE-FG02-91ER40609), the US National Science Foundation, the Bundesministerium fur Bildung und Forschung, Germany, the Alexander von Humboldt Foundation, the UK Engineering and Physical Sciences Research Council, the Polish State Committee for Scientific Research (2 P03B 13820, 02615, 01716, and 09916), the Hungarian Scientific Research Foundation (T14920 and T23790), the EC Marie Curie Foundation, and the Polish-German Foundation.

1