text:
# Longitudinal polarization of hyperons in high 𝑝_⟂ jets in singly polarized 𝑝𝑝 collisions at high energies

## I Introduction

Spin effects in fragmentation processes play an important role in understanding the spin physics in strong interaction. With the developments in experiments, much attention has been attracted to the studies of the different aspects in this connection (see, e.g. Ref.\[1-References\] and references therein). Data are now available from different lepton-induced reaction experiments \[1-References\] and many of the theoretical studies \[References-References\] also concentrate on those reactions. It should be emphasized that the spin effects in the fragmentation processes can also be studied in polarized hadron-hadron collisions by measuring the properties of the hadrons in high transverse momentum $`p_{}`$ jets created in these reactions \[References-References\]. Here, because of the validity of factorization theorem, the hadrons can be regarded as the pure results from the fragmentation of the scattered quark whose polarization can be determined using perturbation theory in quantum chromodynamics (QCD) and the quark distribution functions which determine the polarization before the hard scattering.

Compared with those in the lepton-induced reactions, the study of hadrons in high $`p_{}`$ jets in polarized hadron-hadron collisions has the following properties. (1) Since the involved hard scatterings are strongly interacting processes rather than the electroweak processes involved in the lepton-induced reactions, the corresponding cross section should in general be larger. Furthermore, the luminosity of the incoming proton beams can in general be made higher than that for leptons in particular much higher than that for neutrino. Hence, the statistics can be improved in particular compared with the lepton-induced reactions with neutrino beams. (2) There are many different hard subprocesses which contribute to the production of the high $`p_{}`$ jets. In particular, we have a large contribution from gluon interactions. This, on the one hand, makes the study more complicated than those in lepton-induced reactions, on the other, makes it more interesting since we can use it to study not only quark but also gluon polarization and fragmentation which are even poorly studied yet. Hence, it is also very important to know whether there exist kinematic regions where gluon fragmentation is important or where quark fragmentation dominates. In this case, we can study them separately.

In this paper, we first make an analysis of the contributions of the different hard subprocesses to the production of different hyperons in different kinematic regions in $`pp`$ collisions using a Monte Carlo event generator. Our results show that, there do exist kinematic regions where contributions from quark fragmentation dominate, while in the other gluon fragmentation plays an important role. We then calculate the longitudinal polarization of the hyperons in the regions where quark fragmentation dominates using the same method as we used in studying quark fragmentation in the lepton-induced reactions. In section 2, we present the results for the Monte Carlo analysis of the contributions of different hard processes. In section 3, we briefly summarize the calculations and the main results for the polarizations of the scattered partons from different hard subprocesses. In section 4, we present the results for the polarizations of different hyperons from high $`p_{}`$ jets. A brief summary of the results is given in section 5.

## II Contributions of different hard subprocesses to hyperon production in high $`p_{}`$ jets

We consider the inclusive production of hyperons with high transverse momentum $`p_{}`$ in $`pp`$ collisions at high energy. Here, the factorization theorem can be applied, hence the cross section can be written as,

$$E\frac{d\sigma }{d^3p}(ppHX)=\underset{abcd}{}𝑑x_a𝑑x_bf_a(x_a,\mu ^2)f_b(x_b,\mu ^2)\frac{1}{\pi z_c}\frac{d\widehat{\sigma }}{d\widehat{t}}(abcd)D_c^H(z_c,\mu ^2).$$
(1)
Here, $`f_a(x_a,\mu ^2)`$ and $`f_b(x_b,\mu ^2)`$ are the unpolarized distribution functions of partons $`a`$ and $`b`$ in proton at the scale $`\mu `$, $`x_a`$ and $`x_b`$ are the corresponding momentum fractions carried by $`a`$ and $`b`$; $`D_c^H(z_c,\mu ^2)`$ is the fragmentation function of parton $`c`$ into hadron $`H`$, $`z_c`$ is the momentum fraction of parton $`c`$ carried by the produced $`H`$; $`d\widehat{\sigma }/d\widehat{t}`$ is the cross section at the parton level and $`\widehat{t}`$ is the usual Mandelstam variable of the parton scattering. The cross section of the hard subprocess can be calculated using perturbative QCD. The summation in Eq. (1) runs over all possible subprocesses. To the leading order in pQCD, there are following subprocesses that contribute here, i.e. the scattering between quark(s) and/or antiquark(s), $`qqqq`$, $`q\overline{q}q\overline{q}`$, $`\overline{q}\overline{q}\overline{q}\overline{q}`$; the scattering between gluons or that between gluon and quark, $`qgqg`$, $`\overline{q}g\overline{q}g`$, $`gggg`$; the annihilation processes, $`q\overline{q}gg`$, $`ggq\overline{q}`$. The scattering matrix elements for these subprocesses can easily be calculated and the results can e.g. be found in Ref. \[References\]. For scattering angles less than $`\pi /2`$ in the center of mass frame of the scattering partons, the cross section for $`gggg`$ is the largest, and $`gqgq`$ is the second largest one. They are much higher than the $`qq`$ or $`q\overline{q}`$ scattering in particular at the scattering near $`\pi /2`$. Since the energies of gluons in protons are in general smaller compared with those for the quarks and/or antiquarks, we expect that $`gg`$ scattering plays an important role for jets with moderately large transverse momenta $`p_{}`$. But, for very high $`p_{}`$, $`qq`$ scattering should dominate.

Taking the fragmentation into account, we can estimate the contribution from each subprocess to the production of different hyperons. Presently, this can be conveniently done using Monte Carlo event generators which can give a good fit to the unpolarized data. We used pythia in our calculations and obtained the results for $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^+`$ as a function of pseudorapidity $`\eta `$ in Figs.1 and 2. From these figures, we explicitly see that, at high $`p_{}`$, e.g., $`p_{}`$$`>`$8 GeV at $`\sqrt{s}`$=200 GeV or $`p_{}`$$`>`$13 GeV at $`\sqrt{s}`$=500 GeV, only a few subprocesses dominate. The contributions of quark scatterings $`qqqq`$ (which denotes all the subprocesses in which the initial and final partons are all quarks) and $`qgqg`$ are important whereas that from $`gggg`$ is suppressed significantly. We also calculated the contributions from quark jets or gluon jets. The results are shown in Figs.3 and 4. We see that the contribution from gluon fragmentation is negligibly small in this kinematic region, especially in the large rapidity region (e.g. $`\eta `$$`>`$1.5). Therefore, the uncertainty of the calculations of hyperon polarizations in such kinematic regions caused by the unawareness of the polarized gluon fragmentation function should be very small. In our calculations which will be presented in Section 4, we will just neglect this contribution. However, at lower $`p_{}`$, the contribution from $`gggg`$ or that from gluon fragmentation is very important. As an example, we show in Fig. 5 the corresponding results for $`p_{}`$$`>`$3 GeV at $`\sqrt{s}`$=500 GeV. We see that the contribution of gluon fragmentation is higher than that of quark in the central rapidity regions. Hence, we can study the spin transfer in quark fragmentation in the high $`p_{}`$ region and that in gluon fragmentation in relatively small $`p_{}`$ region.

## III Spin transfer and polarizations of the scattered quarks in different subprocesses

Since we are considering only the elementary subprocess with high $`p_{}`$, the scattering matrix elements can be calculated to high precision using pQCD. The spin transfer factor for each subprocess can be obtained from the helicity amplitudes. The results show that, for all the different kinds of processes $`abcd`$ (Here, $`a`$, $`b`$, $`c`$ and $`d`$ denote different partons, where $`a`$ denotes the one from the longitudinally polarized incoming proton), the longitudinal polarization can be transferred to the outgoing parton $`c`$ and $`d`$. The polarization transfer factors are in general different for different scattering processes. But they have a common feature that they are only functions of a single variable $`y`$, which is defined as $`yk_b(k_ak_c)/k_ak_b`$, where $`k_a`$, $`k_b`$, $`k_c`$ and $`k_d`$ are the four-momenta of the partons $`a`$, $`b`$, $`c`$ and $`d`$ respectively. We summarize these results in Table 1.

From the table, we see the following: (i) Since we neglect quark mass, helicity conservation is valid. This can be seen from the results that $`D_{(3)}^{ac}=D_{(5)}^{ac}=1`$, and that $`D_{(4)}^{ac}=D_{(4)}^{ad}`$. (ii) In the scattering processes $`qq^{}qq^{}`$ and $`qgqg`$, the polarization transfer from one of the colliding objects to the other is the same, and they are also the same as that in the QED process $`e^{}\mu ^{}e^{}\mu ^{}`$ or $`e^{}qe^{}q`$. The result is a monotonously increasing function of $`y`$, which increases from zero to one as $`y`$ goes zero to one. To get a feeling of all these spin transfer factors, we show them as functions of $`y`$ in Fig. 6.

Multiplied by the polarizations of the incoming quarks or gluons, which are determined by the helicity distribution functions of the partons in proton, we obtain the polarizations of the outgoing quarks or gluons. They can be used to calculate the hyperon polarization produced in the fragmentation of these partons. Clearly, the obtained results depend on the polarized parton distribution functions we use. We will also show how strongly the final hyperon polarizations depend on the different sets of parameterization  of these distributions.

## IV Longitudinal hyperon polarization in $`ppHX`$ in high $`p_{}`$ regions

As we have seen in Sec. 2, for hyperons with considerably high transverse momenta, quark fragmentation dominates whereas gluon fragmentation can be neglected. In such kinematic regions, we can calculate the hyperon polarization using the same method as we used in treating $`q_f^0H_i+X`$ in lepton-induced reactions. Measurements of the hyperon polarizations in these kinematic regions can be used as further tests of the different pictures for spin transfer in high energy quark fragmentation processes. The experimental studies can also be made in other kinematic regions where gluon fragmentation plays an important role. By comparing the obtained results with those obtained in the regions where quark fragmentation dominates, we can get useful information for gluon fragmentation. In this section, we present the calculations and the main results in the kinematic regions where quark fragmentation dominates.

### IV.1 The calculation method

The method of calculating the longitudinal polarization $`P_{H_i}`$ of different hyperon $`H_i`$ in the fragmentation of a longitudinally polarized quark $`q_f^0`$ was outlined in Ref. \[References\] using the inclusive process $`e^+e^{}H_i+X`$ as an example. It has been applied to $`e^+e^{}`$ annihilations and deeply inelastic scattering using different lepton beams. We now summarize the main points in the following not only for completeness but also for the following reasons: Here, we would like to emphasize what kinds of inputs that we need in the calculations and how we obtain these inputs. In this way, we can see from where the theoretical uncertainties originate in the calculations and how we should reduce these uncertainties by doing measurements for different hyperons and/or in different kinematic regions.

We consider $`q_f^0H_i+X`$ and divide the produced $`H_i`$’s into the following groups: (a) directly produced and contain the $`q_f^0`$’s; (b) decay products of heavier hyperons which were polarized before their decays; (c) directly produced but do not contain the $`q_f^0`$’s; (d) decay products of heavier hyperons which were unpolarized before their decays. Obviously, hyperons from (a) and (b) can be polarized while those from (c) and (d) are not. Hence, the polarization of $`H_i`$ is given by,

$$P_{H_i}=\frac{\underset{f}{}t_{H_i,f}^FP_f^{(q)}n_{H_i,f}^a+\underset{j}{}t_{H_i,H_j}^DP_{H_j}n_{H_i,H_j}^b}{n_{H_i}^a+n_{H_i}^b+n_{H_i}^c+n_{H_i}^d}.$$
(2)
The different quantities here are defined and obtained in the following way:

(i) $`P_f^{(q)}`$ is the polarization of $`q_f^0`$ which is determined by the initial conditions and the elementary hard scattering processes.

(ii) $`n_{H_i,f}^a`$ is the average number of $`H_i`$’s which are directly produced and contain $`q_f^0`$ of flavor $`f`$, and $`n_{H_i,H_j}^b`$ is that from the decay of $`H_j`$’s which were polarized; $`P_{H_j}`$ is the polarization of $`H_j`$; $`n_{H_i}^a`$, $`n_{H_i}^b`$, $`n_{H_i}^c`$ and $`n_{H_i}^d`$ are the average numbers of $`H_i`$’s in group (a), (b), (c) and (d) respectively. These average numbers of the hyperons of different origins are determined by the hadronization mechanisms and should be independent of the polarization of the initial quarks. Hence, we can calculate them using a hadronization model which give a good description of the unpolarized data. We used Lund model implemented by jetset or lepto or pythia in our calculations for different reactions.

(iii) $`t_{H_i,f}^F`$ is the probability for the polarization of $`q_f^0`$ to be transferred to $`H_i`$ in group (a), i.e. to the $`H_i`$ which contains $`q_f^0`$, and is called the polarization transfer factor, where the superscript $`F`$ stands for fragmentation. It is taken as the fraction of spin carried by the $`f`$-flavor-quark divided by the average number of valence-quark of flavor $`f`$ in $`H_i`$. This is the place where different pictures come in, since the contributions to the hyperon spin from different flavors are different in the SU(6) or the DIS picture. In the SU(6) picture, these contributions can be obtained from the SU(6) wave functions of the hyperons. In the DIS picture, for the $`J^P=(1/2)^+`$ octet hyperons, they are obtained from the DIS data on the spin dependent structure functions and those on hyperon decay. But, for the decuplet hyperons, they are unknown yet. In our calculations, we make a rough estimations of these contributions by taking them as the same as those in the SU(6) picture. It will be interesting and important to study theoretically which picture is applicable in determing $`t_{H_i,f}^F`$. Discussions along this line can be found, for example, in \[References\] and the references given there.

(iv) $`t_{H_i,H_j}^D`$ is the probability for the polarization of $`H_j`$ to be transferred to $`H_i`$ in the decay process $`H_jH_i+X`$ and is called decay polarization transfer factor, where the superscript $`D`$ stands for decay. It is determined by the decay process and is independent of the process in which $`H_j`$ is produced. For the octet hyperon decays, they are extracted from the materials in Review of Particle Properties. But, for the decuplet hyperons, we have to use an estimation based on the SU(6) quark model. (The details are given in e.g \[References,References\]). Hence, we see that decuplet hyperon decay is a major source of the theoretical uncertainties in our calculations of the final $`P_{H_i}`$’s in different reactions. For them, in the DIS picture, neither $`t_{H_i,f}^F`$ nor $`t_{H_i,H_j}^D`$ is known yet. To reduce this uncertainty, it is important to consider the hyperons to which decay contributions are small.

We note also that $`t_{H_i,f}^F`$ and $`t_{H_i,H_j}^D`$ should be considered as universal in the sense that they are the same in different reactions. The differences between different processes come from the polarization $`P_f^{(q)}`$ of the fragmenting quark $`q_f^0`$ and the average numbers $`n_{H_i,f}^a`$, $`n_{H_i,H_j}^b`$, $`n_{H_i}^c`$, and $`n_{H_i}^d`$ for hyperons $`H_i`$ from the different origins. In $`pp`$ collisions that we are considering in this paper, $`P_f^{(q)}`$ is determined by the polarization of the partons before the hard scattering and spin transfer factors in the scattering which are given and discussed in last section.

The average numbers $`n_{H_i,f}^a`$, $`n_{H_i,H_j}^b`$, $`n_{H_i}^c`$, and $`n_{H_i}^d`$ are determined by the hadronization mechanism and the relative abundance of different flavors that take part in the reactions. They should be independent of the process in which $`q_f^0`$ is produced and they are also independent of the polarization of $`q_f^0`$. In $`pp`$ collisions, as we can see from Eq.(1), these average numbers are determined by three factors, i.e., the parton distribution functions $`f_a(x_a,\mu ^2)`$ and $`f_b(x_b,\mu ^2)`$, the cross sections of hard subprocesses $`d\widehat{\sigma }/d\widehat{t}`$, and the fragmentation functions $`D_c^H(z_c,\mu ^2)`$. The cross section $`d\widehat{\sigma }/d\widehat{t}`$ can be calculated with high accuracy from pQCD. For $`f_a(x_a,\mu ^2)`$, $`f_b(x_b,\mu ^2)`$, and $`D_c^H(z_c,\mu ^2)`$, although neither of them is theoretically clear yet, the form of the unpolarized parton densities in nucleon and that of the unpolarized fragmentation functions are empirically known to reasonably high accuracy. We can calculate them using the parameterizations of the structure functions and the phenomenological hadronization models. Since we are not concerned with the different correlations among the produced particles, different hadronization models lead essentially to the same results. Presently, these results can be obtained conveniently from the event generators using Monte Carlo method. We use the Lund string fragmentation model implemented by pythia in our calculations for $`pp`$ collisions.

### IV.2 $`\mathrm{\Lambda }`$ polarization in $`pp\mathrm{\Lambda }X`$ with longitudinally polarized beam

Most of the studies now available on the spin transfer concentrate on $`\mathrm{\Lambda }`$ polarization since $`\mathrm{\Lambda }`$ is most abundantly produced among all the $`J^P=(1/2)^+`$ hyperons and it’s polarization can be measured easily through the decay channel $`\mathrm{\Lambda }p\pi ^{}`$. Meanwhile, as we have seen already in $`e^+e^{}`$ annihilation  and deeply inelastic lepton-nucleon scattering , the origins of $`\mathrm{\Lambda }`$ are also the most complicated. They can be not only directly produced from the fragmentation of $`u`$, $`d`$, or $`s`$ quark, but also from the decays of many different heavier hyperons, such as, $`\mathrm{\Sigma }^0`$, $`\mathrm{\Xi }^{0,}`$, $`\mathrm{\Sigma }^{}(1385)`$, and $`\mathrm{\Xi }^{}(1530)`$. In the large $`p_{}`$ regions in $`pp`$ collisions, the fragmentations of the quarks from the subprocesses $`qgqg`$ and $`qqqq`$ dominate. In particular for hyperons with large rapidities, the scatterings of quarks with large fractional momenta play the most important roles. Since at large $`x`$ in proton, $`u`$ quark density is much higher than that of $`d`$, and both of them are much higher than that of $`s`$, we expect that $`u`$ and $`d`$, particularly $`u`$, quark fragmentation dominates the $`\mathrm{\Lambda }`$ production in the high $`p_{}`$ and large rapidity region. However, $`u`$ or $`d`$ quark carries no spin of $`\mathrm{\Lambda }`$ in the SU(6) picture and only a small fraction of the spin of $`\mathrm{\Lambda }`$ in the DIS picture, we expect that, just as in deeply inelastic scattering , the magnitudes of the $`P_\mathrm{\Lambda }`$ obtained in both models should be quite small, much smaller than those obtained in $`e^+e^{}`$ annihilation near the $`Z^0`$ pole. The influence of heavier hyperon decay should be very important and the induced theoretical uncertainties in the calculations should be relatively large.

Using the event generator pythia, we obtained the average numbers $`n_{\mathrm{\Lambda },f}^a`$, $`n_{\mathrm{\Lambda },H_j}^b`$, $`n_\mathrm{\Lambda }^c`$, and $`n_\mathrm{\Lambda }^d`$ for $`\mathrm{\Lambda }`$ production in $`pp\mathrm{\Lambda }X`$. By inserting them into Eq. (2), we obtain $`P_\mathrm{\Lambda }`$ as a function of the pseudorapidity $`\eta `$. We show the different contributions to $`\mathrm{\Lambda }`$ production at $`\sqrt{s}`$ =500 GeV and $`p_{}`$$`>`$13 GeV as functions of $`\eta `$ in Fig. 7 and $`P_\mathrm{\Lambda }`$ in Fig. 8. The moving direction of polarized beam proton is taken as the rapidity axis. From these results, we see in particular the following characteristics.

First, from Fig.7, we see that the contributions from $`u`$-quark fragmentation indeed dominate the $`\mathrm{\Lambda }`$’s, in particular at large $`\eta `$. There are several direct consequences following from this : (i) the contributions from heavier hyperon decays are large, even larger than the directly produced contributions. To show this point more clearly, we present as an example in Fig. 9 the different contributions to $`P_\mathrm{\Lambda }`$ in the SU(6) picture. We see that the magnitudes of the contributions to $`P_\mathrm{\Lambda }`$ from $`\mathrm{\Sigma }^0`$ and $`\mathrm{\Sigma }^0`$ decays are both much larger than those from directly produced $`\mathrm{\Lambda }`$’s. We see that even the qualitative tendency are influenced very much by these decay contributions. (ii) the resulting $`|P_\mathrm{\Lambda }|`$ is relatively small, it is less than 3%. (iii) the differences between the results from SU(6) and those from the DIS picture are small.

Second, the contributions from those hyperons which are directly produced and contain the fragmenting $`q_f^0`$’s \[origin (a) mentioned in section 4.A\] or from the decays of the heavier hyperons which are directly produced and contain the fragmenting $`q_f^0`$’s \[the type of heavier hyperons of origin (a)\] are much larger than the corresponding contributions in the lepton-induced reactions. This is very interesting since these hyperons are all polarized and the polarizations are proportional to that of $`q_f^0`$. In fact, in the limiting case that there are only these two types of hyperons contribute, we have,

$$P_{H_i}^{(lim)}=\frac{\underset{f}{}t_{H_i,f}^FP_f^{(q)}n_{H_i,f}^a+\underset{j}{}t_{H_i,H_j}^DP_{H_j}^{(a)}n_{H_i,H_j}^b}{n_{H_i}^a+n_{H_i}^b}.$$
(3)
Here, we use the superscript $`(a)`$ in $`P_{H_j}^{(a)}`$ to denote the polarization of hyperon $`H_j`$ of only origin (a). It is clear that, from Eq.(2), for such heavier hyperons,

$$P_{H_j}^{(a)}=\frac{\underset{f}{}t_{H_j,f}^FP_f^{(q)}n_{H_j,f}^a}{n_{H_j}^a}.$$
(4)
By inserting it into Eq.(3), we obtain that,

$$P_{H_i}^{(lim)}=\frac{\underset{f}{}[t_{H_i,f}^Fn_{H_i,f}^a+\underset{j}{}t_{H_i,H_j}^Dt_{H_j,f}^Fn_{H_j,f}^aR_{H_i,H_j}]P_f^{(q)}}{n_{H_i}^a+\underset{j}{}n_{H_j}^aR_{H_i,H_j}},$$
(5)
where $`R_{H_i,H_j}`$ denotes the branching ratio for the decay $`H_jH_i+X`$. We see clearly that, in this limiting case, the $`P_{H_i}`$ reflects more directly the polarization of the fragmenting quark $`q_f^0`$. It is proportional to the $`P_f^{(q)}`$’s and the proportional constants depend on the picture for spin transfer and also the decay polarization transfer factors involved. We see clearly that $`P_\mathrm{\Lambda }`$ should depends rather strongly on $`P_f^{(q)}`$, which is determined by the spin transfer factors given in Table 1 and the polarized parton distribution functions. To see how sensitive it is, we used two different sets of parameterizations in the calculations. We see that the resulting difference in $`P_\mathrm{\Lambda }`$ is indeed quite significant.

### IV.3 Polarization of other hyperons

From the results that we obtained in the last subsection, we see that the contributions from the decay of heavier hyperons to $`\mathrm{\Lambda }`$ are large and the $`\mathrm{\Lambda }`$ polarization is very small. Since the contributions of different flavors to the spins of the $`J^P=(3/2)^+`$ hyperons are unknown yet in the DIS picture and only models are known to calculate the spin transfer in their decay processes, there are considerably large theoretical uncertainties in the results that we can obtain yet. On the other hand, as we have already seen in $`e^+e^{}`$$``$$`Z^0`$$``$$`HX`$ and in $`\mu ^{}p\mu ^{}HX`$, the origins of other $`J^P=(1/2)^+`$ hyperons, i.e., $`\mathrm{\Sigma }`$ and $`\mathrm{\Xi }`$, are rather simple. The decay contributions are much smaller than the corresponding contributions for $`\mathrm{\Lambda }`$ production. The theoretical uncertainties in this case can be significantly reduced. In particular in the large rapidity region, the $`u`$ quark fragmentation plays the dominant role and $`u`$ carried most of the hyperon $`\mathrm{\Sigma }^+`$’s spin. Hence, the resulting $`\mathrm{\Sigma }^+`$ polarization in $`pp\mathrm{\Sigma }^+X`$ should be much larger than that for $`\mathrm{\Lambda }`$. Furthermore, also because of the dominance of the $`u`$ quark fragmentation, the production rate of $`\mathrm{\Sigma }^+`$ should be comparable with that of $`\mathrm{\Lambda }`$, which implies that the statistics for studying $`\mathrm{\Sigma }^+`$ should be similar to that of $`\mathrm{\Lambda }`$.

Using the event generator pythia, we calculate the different contributions to the inclusive $`\mathrm{\Sigma }^+`$ production in $`pp`$ collisions. Using Eq.(2), we obtain $`P_{\mathrm{\Sigma }^+}`$ as a function of $`\eta `$ at $`p_{}`$$`>`$13 GeV. The results are shown in Figs. 10 and 11 (a). Here, we see that the decay contribution from heavier hyperons is very small. It takes only a few percents for $`p_{}`$$`>13`$ GeV at $`\sqrt{s}`$=500 GeV. We see also that the obtained $`P_{\mathrm{\Sigma }^+}`$ is much larger than $`P_\mathrm{\Lambda }`$ and the difference between those results based on the two different pictures are much more significant, in particular for large $`\eta `$. From these calculations, we also see again that the contributions from the directly produced $`\mathrm{\Sigma }^+`$’s which contain the fragmenting quark play the dominate role. In fact, they give more than 75% of the whole produced $`\mathrm{\Sigma }^+`$ at $`p_{}`$$`>`$$`13`$ GeV. In particular for $`\eta `$$`>`$$`1.5`$, those directly produced and contain the fragmenting $`u`$-quarks take more than 93%. If we neglect all other contributions, we obtain the $`P_{\mathrm{\Sigma }^+}`$ in this limiting case, i.e. the $`P_{\mathrm{\Sigma }^+}`$ for those only from the origin (a) as,

$$P_{\mathrm{\Sigma }^+}^{(lim)}=P_{\mathrm{\Sigma }^+}^{(a)}=t_{\mathrm{\Sigma }^+,u}^FP_u^{(q)}.$$
(6)
We see that it is directly proportional to the polarization of the $`u`$-quark, $`P_u^{(q)}`$, after the hard scattering and the proportional constant is just the fragmentation spin transfer $`t_{\mathrm{\Sigma }^+,u}^F`$. This provides indeed an excellent place to test different pictures for spin transfer in fragmentation processes. It may provide also a nice tool to study the polarization of the quarks. Further studies along this line are under way. From Fig.11, we see also that the results for $`P_{\mathrm{\Sigma }^+}`$ from two different sets of parameterizations of the helicity parton distributions are rather small. This is because $`P_{\mathrm{\Sigma }^+}`$ depends mainly on $`\mathrm{\Delta }u`$, which is essentially the same in the two sets of parameterizations.

We also make similar calculations for other hyperons, i.e., $`\mathrm{\Sigma }^{}`$, $`\mathrm{\Xi }^0`$ and $`\mathrm{\Xi }^{}`$. The results are shown in Fig. 11 (b)-(d). We see that, $`|P_\mathrm{\Sigma }^{}|`$ exhibits similar behavior as $`|P_{\mathrm{\Sigma }^+}|`$ but has opposite sign. This is because, similar to $`\mathrm{\Sigma }^+`$ where $`u`$-quark fragmentation dominates, $`d`$-quark fragmentation dominates $`\mathrm{\Sigma }^{}`$ production. But $`\mathrm{\Delta }d`$ has opposite sign as $`\mathrm{\Delta }u`$ in proton. Both $`P_{\mathrm{\Xi }^0}`$ and $`P_\mathrm{\Xi }^{}`$ are determined mainly by the fragmentation of $`s`$-quark but with some influence from $`u`$ or $`d`$. Since $`u`$ or $`d`$ contributes negatively to the spin of $`\mathrm{\Xi }`$, and $`\mathrm{\Delta }u`$ is positive in proton but $`\mathrm{\Delta }d`$ is negative, we should have a negative of $`u`$ quark fragmentation to $`P_{\mathrm{\Xi }^0}`$ but a positive contribution from $`d`$ to $`P_\mathrm{\Xi }^{}`$. This makes $`P_{\mathrm{\Xi }^0}`$ and $`P_\mathrm{\Xi }^{}`$ behave a bit different. The difference can be seen in Fig.11. All these different features can be used to make further tests to the SU(6) or DIS picture used for the spin transfer in fragmentation processes.

## V Summary

In summary, after a Monte Carlo analysis of the different contributions, we calculated the longitudinal polarizations for different hyperons in high $`p_{}`$ jets in $`pp`$ collisions with longitudinally polarized beams using two different pictures: the SU(6) and the DIS pictures. The results show that the $`\mathrm{\Lambda }`$ polarization is rather small and has large contributions of different heavier hyperons decay. On the contrary, the contribution to $`\mathrm{\Sigma }^+`$ production or that to $`\mathrm{\Sigma }^{}`$ production is rather pure and the polarization is much larger. Hence, the measurements of these polarizations can give a good test to the validity of the two pictures. Such measurements can e.g. be carried out at RHIC in the near future.

## Acknowledgments

We would like to thank Li Shi-yuan and Xie Qu-bing and other members in the theoretical particle physics group of Shandong University for helpful discussions. This work was supported in part by the National Science Foundation of China (NSFC) and the Education Ministry of China under Huo Ying Dong Foundation.

Table caption

Table 1: Polarization transfer factors $`D_{(i)}^{ac}(y)`$ and $`D_{(i)}^{ad}(y)`$ in the $`i`$-th kind of the elementary processes $`abcd`$ when the incoming parton $`a`$ is longitudinally polarized. Here, $`A`$ and $`B`$ are defined as $`A310y+13y^26y^3+3y^4`$ and $`B32y+y^22y^3+3y^4`$.

Figure captions

Fig.1: Contributions of the different parton level subprocesses to $`\mathrm{\Lambda }`$ production in $`pp\mathrm{\Lambda }X`$.

Fig.2: Contributions of the different parton level subprocesses to $`\mathrm{\Sigma }^+`$ production in $`pp\mathrm{\Sigma }^+X`$.

Fig.3: The contributions of quark and gluon jets to $`\mathrm{\Lambda }`$ production in $`pp\mathrm{\Lambda }X`$.

Fig.4: The contributions of quark and gluon jets to $`\mathrm{\Sigma }^+`$ production in $`pp\mathrm{\Sigma }^+X`$.

Fig.5: The contributions of quark and gluon jets to $`\mathrm{\Lambda }`$ production (a) and $`\mathrm{\Sigma }^+`$ production (b) at $`p_{}`$$`>`$3 GeV in $`ppHX`$.

Fig.6: Polarization transfer factors $`D_{(i)}^{ac}(y)`$ and $`D_{(i)}^{ad}(y)`$ in the $`i`$-th kind of the elementary processes $`abcd`$ in Table 1.

Fig.7: Different contributions to $`\mathrm{\Lambda }`$ production with $`p_{}`$$`>`$13 GeV in events originating from the $`d`$, $`u`$, or $`s`$ quark fragmentation as a function of $`\eta `$ in $`pp\mathrm{\Lambda }X`$ at $`\sqrt{s}=500`$ GeV. In (a)-(c), we see the five types of contributions from the fragmentation of $`d`$, $`u`$, or $`s`$ quark respectively. The solid line denotes those which are directly produced and contain the fragmenting quark; the dashed, dash-dotted, upper dotted, lower dotted, denote those from the decay of $`\mathrm{\Sigma }^0`$, $`\mathrm{\Xi }`$, $`\mathrm{\Sigma }^{}`$ and $`\mathrm{\Xi }^{}`$ which contain the fragmenting quark. In (d)-(f), we show the fractions $`R_i`$ which are the ratios of the corresponding contributions to the sums of all different contributions.

Fig.8: Longitudinal $`\mathrm{\Lambda }`$ polarization $`P_\mathrm{\Lambda }`$ as a function of $`\eta `$ for $`\mathrm{\Lambda }`$ with $`p_{}`$$`>`$13 GeV in $`pp\mathrm{\Lambda }X`$at $`\sqrt{s}=`$ 500 GeV with one of the beams is longitudinally polarized. ‘GRSV2000/GRV98’ denotes the Standard LO set of GRSV2000 for polarized distribution functions and GRV98 LO set for unpolarized distribution functions we used here; ‘GRSV96/GRV94’ denotes the Standard LO set of GRSV96 and GRV94 LO set correspondingly.

Fig.9: Different contributions to $`P_\mathrm{\Lambda }`$ in the SU(6) picture for $`\mathrm{\Lambda }`$ with $`p_{}`$$`>`$13 GeV in $`pp\mathrm{\Lambda }X`$ at $`\sqrt{s}=`$ 500 GeV in the case that the Standard LO set of GRSV2000 for polarized distribution functions and GRV98 LO set for unpolarized distribution functions are used.

Fig.10: Different contribution to $`\mathrm{\Sigma }^+`$ production with $`p_{}`$$`>`$13 GeV in $`pp\mathrm{\Sigma }^+X`$ at $`\sqrt{s}=`$ 500 GeV. The solid and dashed lines are respectively the contributions which are directly produced and contain the fragmenting quark $`u`$ and $`s`$ respectively; the dotted line corresponds to those which are originate from the decay of polarized heavier hyperons.

Fig.11: Polarizations for $`\mathrm{\Sigma }^+`$, $`\mathrm{\Sigma }^{}`$, $`\mathrm{\Xi }^0`$ and $`\mathrm{\Xi }^{}`$ with $`p_{}`$$`>`$13 GeV in polarized $`ppHX`$ as a function of $`\eta `$ at $`\sqrt{s}=`$ 500 GeV. Two sets of polarized distribution functions as same as for $`\mathrm{\Lambda }`$ are used.


format: text

text:
# Master equation approach to the conjugate pairing rule of Lyapunov spectra for many-particle thermostatted systems

## I Introduction

The Lyapunov exponent is an essential concept to express the instability of orbits in a dynamical system. It is introduced as an exponential expansion (or contraction) rate of an infinitesimal perturbation of orbits, and its positivity implies that the system is chaotic. In general there is a Lyapunov exponent for each independent direction of the infinitesimal perturbation of the orbit, and the sorted set of such Lyapunov exponents is called the Lyapunov spectrum, and has been the subject of study in many-particle systems. For example, the existence of its thermodynamic limit , an effect of the rotational degrees of freedom of molecules , its stepwise structure and the Lyapunov modes , and a tracer particle effect  have been observed and discussed in the Lyapunov spectra of many-particle chaotic systems.

A characteristic of the Lyapunov spectrum that is known in Hamiltonian systems, is that the Lyapunov exponents appear as a pair, namely any positive Lyapunov exponent accompanies a negative Lyapunov exponent with its opposite sign . This characteristic, which is based on the symplectic structure of the Hamiltonian mechanics, is not correct in non-Hamiltonian systems, but it is interesting to know how it is modified in quasi-Hamiltonian systems such as a Hamiltonian system coupled to a thermodynamic reservoir. This problem has been considered in some thermostatted dynamics with a term to extract the heat produced in the system by external force fields, and led to the proposal of the conjugate pairing rule for thermostatted systems, which claims that the sum of any Lyapunov exponent pair is not zero but a constant regardless of the exponent numbers . This conjecture was confirmed by many following numerical calculations . This also led to the discovery that some thermostatted systems contain hidden Hamiltonian structure . The pairing rule is not only interesting as a mathematical structure of the thermostatted system but also valuable for a practical use; the conjugate pairing rule for the thermostatted system allows us to calculate non-equilibrium transport coefficients (eg. conductivity and viscosity) from only one pair of the Lyapunov exponents, such as the largest and smallest Lyapunov exponents only .

A problem is that the necessary and the sufficient conditions for the conjugate pairing rule to hold for thermostatted systems is not clearly known. The conjugate pairing rule for the iso-kinetic thermostatted system with a color field was proved for the soft core interaction potential  and the hard core interaction potential , regardless of the number of particles. A similar discussion was done in Nosé Hamiltonian dynamics . These works give the sufficient conditions for the conjugate pairing rule. On the other hand, it was suggested numerically that it can be violated in the presence of a magnetic field  and in inhomogenously thermostatted systems such as a system under temperature gradient  or a system in which the peculiar momenta are thermostatted . Another numerical work also suggested that it is not exact in the iso-energetic thermostat with a finite number of particles , although the iso-energetic thermostat should be equivalent to the iso-kinetic thermostat as the number of the particles goes to infinity, namely in the thermodynamic limit . A special interest is the iso-kinetic thermostatted system with a shear field, which is described by the Sllod equation for the planar Coutte flow . Early investigations supported the conjugate pairing rule for such a system . Refs.  suggested a small deviation from the conjugate pairing rule. An analytical consideration showed that the deviation from the conjugate pairing rule should be at most fourth order in the shear rate in the case of a small shear rate in the thermodynamic limit . However a recent numerical calculation with a more careful numbering of the Lyapunov exponents and with numerical error bars showed that within the numerical precision the conjugate pairing rule was satisfied . After all these trials, a justification of the conjugate pairing rule for the iso-kinetic thermostatted system with a shear field still remains as an open problem.

Analytical calculation of the full Lyapunov spectra for many-particle systems is still not an easy matter, and so far full Lyapunov spectra have been calculated mainly using numerical approaches. On the other hand, recently an analytical approach to the full Lyapunov spectra was proposed for many-particle systems with a random interaction . This approach describes the tangent vector dynamics by a master equation, and allows the calculation of all individual positive Lyapunov exponents through the average of the magnitude of the tangent vector. It was used to explain the stepwise structure of the Lyapunov spectrum for many-particle Hamiltonian systems. However a generalization of this approach to non-equilibrium systems, especially to thermostatted systems, is not known yet.

The purpose of this paper is to generalize the master equation approach to the Lyapunov spectrum to non-equilibrium thermostatted systems, and to discuss the conjugate pairing rule of the Lyapunov spectrum for the iso-kinetic thermostatted system with a shear field. Especially we propose a method to calculate negative Lyapunov exponents using a time-reversed master equation. In this paper we concentrate on the case where the particles interact with a Gaussian white randomness. In this case the master equation that describes the tangent vector dynamics is attributed to the Fokker-Planck equation. We also restrict our consideration in the thermodynamic limit, in which the fluctuations of the friction coefficient can be neglected . Under these conditions we show that the conjugate pairing rule for the thermostatted system with a shear field given by the Sllod equation is satisfied. As a special case we also discuss briefly an explicit form to connect the Lyapunov exponents with the time-correlation of the interaction matrix in a thermostatted system without a shear field.

## II Isokinetic Thermostatted System with a Shear Field and its Tangent Vector Dynamics

We consider non-equilibrium systems with an iso-kinetic thermostat. Our consideration includes the case where a shear flow is sustained by an external restriction, and for simplicity we consider a two-dimensional system consisting of $`N`$ particles with the same mass $`m`$. We introduce $`𝐪^{(j)}(t)(q_x^{(j)}(t),q_y^{(j)}(t))^T`$ and $`𝐩^{(j)}(t)(p_x^{(j)}(t),p_y^{(j)}(t))^T`$ as the spatial coordinate vector and the momentum vector of the $`j`$-particle, respectively, at time $`t`$ with the transpose operation $`T`$. (Note that all vectors in this paper are introduced as column vectors.) Equations for $`𝐪^{(j)}(t)`$ and $`𝐩^{(j)}(t)`$ are expressed as

$`{\displaystyle \frac{d𝐪^{(j)}(t)}{dt}}={\displaystyle \frac{1}{m}}𝐩^{(j)}(t)+\gamma \mathrm{\Xi }_2𝐪^{(j)}(t)`$ (1)
$`{\displaystyle \frac{d𝐩^{(j)}(t)}{dt}}={\displaystyle \frac{U(t)}{𝐪^{(j)}(t)}}\left[\gamma \mathrm{\Xi }_2+\alpha (t)I_2\right]𝐩^{(j)}(t)`$ (2)
where $`U(t)`$ is the potential energy as a function of $`𝐪^{(j)}(t),j=1,2,\mathrm{},N`$ and $`t`$ only, and we introduce $`\mathrm{\Xi }_{2k}`$ is the $`2k\times 2k`$ matrix defined by

$`\mathrm{\Xi }_{2k}\left(\begin{array}{cc}0_k& I_k\\ 0_k& 0_k\end{array}\right)`$ (5)
with the $`k\times k`$ identical matrix $`I_k`$ and the $`k\times k`$ null matrix $`0_k`$. Here $`\gamma `$ is the shear rate as an external parameter, namely a constant gradient of the $`x`$ component of the local velocity in the $`y`$ direction, and $`\alpha (t)`$ is defined by

$`\alpha (t){\displaystyle \frac{_{j=1}^N𝐩^{(j)}(t)^T\left(\frac{U(t)}{𝐪^{(j)}(t)}+\gamma \mathrm{\Xi }_2𝐩^{(j)}(t)\right)}{_{j=1}^N|𝐩^{(j)}(t)|^2}}`$ (6)
so that the total kinetic energy is constant in time: $`d\left[_{j=1}^N|𝐩^{(j)}(t)|^2/(2m)\right]/dt=0`$. Eqs. (1) and (2) are called the Sllod equation for the planar Coutte flow with the isokinetic thermostat, and gives the model for the system driven by external fields and/or a shear rate with an attached heat reservoir which removes the energy generated inside the system and maintains the temperature of the system constantly in time. As an example described by Eqs. (1) and (2), other than the system with a shear field, we may mention the color field system in which the system consists of many particles with charges of different signs and is driven by an external electric field .

In general, the quantity $`\alpha (t)`$, which is interpreted as the friction coefficient, depends on the coordinates and the momenta of the particles, so is variable in time. However, it is known that the fluctuation of the quantity $`\alpha (t)`$ is small in a system consisting of many particles . (For a justification of this point by the kinetic approach see Ref. , which shows that the quantity $`\alpha (t)`$ fluctuates with the order of $`1/\sqrt{N}`$ around a fixed value.) Based on this fact, in this paper we consider only the system which consists of enough particles so that the friction coefficient $`\alpha (t)`$ in Eq. (2) can be replaced by a fixed constant $`\overline{\alpha }`$.

For a convenience we represent the $`4N`$-dimensional phase space vector $`𝚪(t)`$ as a vector $`(q_x^{(1)}(t),q_x^{(2)}(t),\mathrm{},q_x^{(N)}(t),`$ $`q_y^{(1)}(t),q_y^{(2)}(t),\mathrm{},q_y^{(N)}(t),`$ $`p_x^{(1)}(t),p_x^{(2)}(t),\mathrm{},p_x^{(N)}(t),`$ $`p_y^{(1)}(t),p_y^{(2)}(t),\mathrm{},p_y^{(N)}(t))^T`$. Using this notation and the assumption explained in the previous paragraph we obtain the equation

$`{\displaystyle \frac{d\delta 𝚪(t)}{dt}}=(t)\delta 𝚪(t)`$ (7)
for the tangent vector $`\delta 𝚪(t)`$. Here the matrix $`(t)`$ is given by

$`(t)\left(\begin{array}{cc}\mathrm{\Phi }& I_{2N}/m\\ R(t)& \mathrm{\Psi }\end{array}\right)`$ (10)
with $`2N\times 2N`$ matrices $`\mathrm{\Phi }`$, $`\mathrm{\Psi }`$ and $`R`$ defined by

$`\mathrm{\Phi }\gamma \mathrm{\Xi }_{2N},`$ (11)
$`\mathrm{\Psi }\gamma \mathrm{\Xi }_{2N}\overline{\alpha }I_{2N},`$ (12)
$`R(t){\displaystyle \frac{^2U(t)}{𝐪(t)𝐪(t)}}`$ (13)
where we introduced $`𝐪(t)`$ as a vector $`(q_x^{(1)}(t),q_x^{(2)}(t),\mathrm{},`$ $`q_x^{(N)}(t),`$ $`q_y^{(1)}(t),q_y^{(2)}(t),\mathrm{},q_y^{(N)}(t))^T`$.

## III Random Interactions and Master Equations for the Tangent Vector Dynamics

In this section we introduce a random interaction between the particles, and obtain the two kinds of master equations corresponding to the time-forward tangent vector dynamics and the time-reversed tangent vector dynamics by using Kramers-Moyal expansion technique.

### A Fokker-Planck equation for the forward dynamics of the tangent vector

We consider the case that each particle interacts with the other particles randomly enough so that the matrix $`R(t)(R_{jk}(t))`$ can be regarded as a Gaussian white random matrix satisfying the conditions

$`R_{\mu _1\nu _1}(t_1)R_{\mu _2\nu _2}(t_2)\mathrm{}R_{\mu _{2n1}\nu _{2n1}}(t_{2n1})=0,`$ (14)
$`R_{\mu _1\nu _1}(t_1)R_{\mu _2\nu _2}(t_2)\mathrm{}R_{\mu _{2n}\nu _{2n}}(t_{2n})`$ (15)
$`={\displaystyle \underset{P_d}{}}D_{\mu _{j_1}\nu _{j_1}\mu _{j_2}\nu _{j_2}}D_{\mu _{j_3}\nu _{j_3}\mu _{j_4}\nu _{j_4}}`$ (16)
$`\mathrm{}D_{\mu _{j_{2n1}}\nu _{j_{2n1}}\mu _{j_{2n}}\nu _{j_{2n}}}`$ (17)
$`\times \delta (t_{j_1}t_{j_2})\delta (t_{j_3}t_{j_4})\mathrm{}\delta (t_{j_{2n1}}t_{j_{2n}})`$ (18)
for any integer $`n`$ and a $`4`$-th rank constant tensor $`D_{jkln}`$, where we take the sum over only the permutation $`P_d:`$ $`(1,2,\mathrm{}`$ $`,2n)`$ $``$ $`(j_1,j_2,\mathrm{},j_{2n})`$, and the bracket $`\mathrm{}`$ means the ensemble average over random processes.

Under the randomness conditions (14) and (18) the time-evolutional equation (7) is regarded as a stochastic equation of the Langevin type, and its corresponding master equation for the the probability density $`\rho ^{(+)}(\delta 𝚪,t)`$ for the tangent vector $`\delta \mathrm{\Gamma }`$ at time $`t`$ is given by

$`{\displaystyle \frac{}{t}}\rho ^{(+)}(\delta 𝚪,t)`$ (19)
$`={\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \frac{}{\delta q_\mu }}\left(\mathrm{\Phi }_{\mu \nu }\delta q_\nu +{\displaystyle \frac{\delta _{\mu \nu }}{m}}\delta p_\nu \right)\rho ^{(+)}(\delta 𝚪,t)`$ (20)
$`{\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \frac{}{\delta p_\mu }}\mathrm{\Psi }_{\mu \nu }\delta p_\nu \rho ^{(+)}(\delta 𝚪,t)`$ (21)
$`+{\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \underset{\mu ^{}=1}{\overset{2N}{}}}{\displaystyle \underset{\nu ^{}=1}{\overset{2N}{}}}{\displaystyle \frac{1}{2}}D_{\mu ^{}\mu \nu ^{}\nu }\delta q_\mu \delta q_\nu `$ (22)
$`\times {\displaystyle \frac{^2}{\delta p_\mu ^{}\delta p_\nu ^{}}}\rho ^{(+)}(\delta 𝚪,t)`$ (23)
applying the Kramers-Moyal expansion technique to the dynamics (7). Here $`\delta q_j`$ and $`\delta p_j`$ are the $`j`$-th components of the coordinate part $`\delta 𝐪`$ and the momentum part $`\delta 𝐩`$ in the tangent vector $`\delta \mathrm{\Gamma }=(\delta 𝐪,\delta 𝐩)^T`$, respectively, and $`\mathrm{\Phi }_{\mu \nu }`$ and $`\mathrm{\Psi }_{\mu \nu }`$ are the matrix elements of the matrix $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ defined by Eqs. (11) and (12), respectively. The derivation of Eq. (23) is given in Appendix A. Eq. (23) in the special case of $`\mathrm{\Phi }=0_{2N}`$ and $`\mathrm{\Psi }=0_{2N}`$ have already been used to discuss the stepwise structure of the Lyapunov spectrum for a many-particle Hamiltonian system .

### B Anti-Fokker-Planck equation for the time-reversed dynamics of the tangent vector

As shown in Ref. , in the case of $`\gamma =0`$ and $`\overline{\alpha }=0`$ the Fokker-Planck equation (23) provides the positive Lyapunov exponents as the time averaged exponential rate of the randomness average by the probability density $`\rho ^{(+)}(𝚪,t)`$ in the time evolution of infinitesimal perturbations of the dynamical variables. However this method does not provide directly the negative Lyapunov exponents, because in the stochastic system the randomness average of the distance between the infinitesimal nearby trajectories should not shrink in the infinite time limit. This was not a problem in the Hamiltonian system discussed in Ref. , because in the Hamiltonian system the absolute values of the negative Lyapunov exponents are the same with the positive Lyapunov exponents . However in the thermostatted system discussed in this paper such a simple relation of the negative and the positive Lyapunov exponents can not be expected any more. In order to overcome this problem and to provide the negative Lyapunov exponents using the master equation approach directly, we use the fact that the negative Lyapunov exponents can be regarded as the positive Lyapunov exponents for the time-reversed motion. This fact has been used in some works to calculate the negative Lyapunov exponents for chaotic systems .

In the iso-kinetic thermostatted system with a shear field the time-reversed motion is expressed by the ”time-reversed mapping” : $`tt`$, $`𝐪𝐪`$, $`𝐩𝐩`$ and $`\gamma \gamma `$. The transformation $`\gamma \gamma `$ is justified by the fact that the direction of the shear flow changes to the opposite direction in the time-reversed motion. This justifies the time-reversal transformation $`\overline{\alpha }\overline{\alpha }`$ for the friction coefficient by Eq. (6). The time-reversed mapping leads to the time-reversed transformations

$`tt`$ (24)
$`\delta 𝐪\delta 𝐪`$ (25)
$`\delta 𝐩\delta 𝐩`$ (26)
$`\mathrm{\Phi }\mathrm{\Phi }`$ (27)
$`\mathrm{\Psi }\mathrm{\Psi }`$ (28)
for the tangent vector dynamics, noting the relations (11) and (12) to connect the matrices $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ with the shear rate $`\gamma `$ and the friction coefficient $`\overline{\alpha }`$. It is important to note that the tensor $`D_{jkln}`$ and the tangent dynamical equation (7) themselves are invariant under the transformations (24-28).

Now we introduce the Fokker-Planck equation for the probability density $`\rho ^{()}(\delta 𝚪,t)`$ for the time-reversed tangent vector at time $`t`$ as the transformed equation of the Fokker-Planck equation (23) by the transformations (24-28), namely

$`{\displaystyle \frac{}{t}}\rho ^{()}(\delta 𝚪,t)`$ (29)
$`={\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \frac{}{\delta q_\mu }}\left(\mathrm{\Phi }_{\mu \nu }\delta q_\nu +{\displaystyle \frac{\delta _{\mu \nu }}{m}}\delta p_\nu \right)\rho ^{()}(\delta 𝚪,t)`$ (30)
$`{\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \frac{}{\delta p_\mu }}\mathrm{\Psi }_{\mu \nu }\delta p_\nu \rho ^{()}(\delta 𝚪,t)`$ (31)
$`{\displaystyle \underset{\mu =1}{\overset{2N}{}}}{\displaystyle \underset{\nu =1}{\overset{2N}{}}}{\displaystyle \underset{\mu ^{}=1}{\overset{2N}{}}}{\displaystyle \underset{\nu ^{}=1}{\overset{2N}{}}}{\displaystyle \frac{1}{2}}D_{\mu ^{}\mu \nu ^{}\nu }\delta q_\mu \delta q_\nu `$ (32)
$`\times {\displaystyle \frac{^2}{\delta p_\mu ^{}\delta p_\nu ^{}}}\rho ^{()}(\delta 𝚪,t).`$ (33)
This is simply the equation with the opposite sign of the diffusion term to the forward Fokker-Planck equation (23), and is interpreted as the master equation to describe the time-evolution of the tangent vector whose initial condition is the time-reversed initial condition to the Fokker-Planck equation (23). We call this equation for the probability density $`\rho ^{()}(\delta 𝚪,t)`$ the ”anti-Fokker-Planck equation” in this paper, and calculate the negative Lyapunov exponents as the exponential rate of the randomness average by the probability density $`\rho ^{()}(\delta 𝚪,t)`$ in the time evolution of infinitesimal perturbations of the dynamical variables in the minus infinite time.

It is important to note the difference between the anti-Fokker-Planck equation (33) and the so called ”backward Fokker-Planck equation” . The backward Fokker-Planck equation describes the dynamics before an initial time in which the initial condition is the same as in the forward Fokker-Planck equation to describe the dynamics after the initial time. However, in order to calculate the negative Lyapunov exponents from the time-reversed motion, we must use the different initial condition which has the opposite sign of the momentum to the initial condition in the forward Fokker-Planck equation. In the equilibrium case expressed by $`\gamma =0`$ and $`\overline{\alpha }=0`$ the backward Fokker-Planck equation coincides with the anti-Fokker-Planck equation, but otherwise it can not be used to calculate the negative Lyapunov exponents.

The anti-Fokker-Planck equation is analogous to the ”anti-Lorentz-Boltzmann equation” which was introduced to calculate the negative Lyapunov exponents using the kinetic approach . In this approach the anti-Lorentz-Boltzmann equation is given as the Lorentz-Boltzmann equation where the collision operator has the opposite sign to the ordinary Lorentz-Boltzmann equation in an equilibrium or a non-equilibrium stationary state.

At least, in the equilibrium case expressed by $`\gamma =0`$ and $`\overline{\alpha }=0`$ the anti-Fokker-Planck equation must provide the negative Lyapunov exponents as the opposites of the positive Lyapunov exponents calculated by using the Fokker-Planck equation (23). In the next section we show it actually as the special case of more general results.

## IV Conjugate Pairing Rule for Thermostatted Systems with a Shear

We have to know the time evolution of the amplitude of the tangent vector in order to calculate the Lyapunov exponents. Such a time evolution for the forward movement (the time-reversed movement) is expressed as the dynamics of the diagonal elements of the matrix $`\mathrm{{\rm Y}}^{(+)}(t)`$ (the matrix $`\mathrm{{\rm Y}}^{()}(t)`$) given by

$`\mathrm{{\rm Y}}^{(\pm )}(t)\delta 𝐪\delta 𝐪^T_t^{(\pm )}.`$ (34)
Here the bracket $`\mathrm{}_t^{(\pm )}`$ means to take the average by the probability density $`\rho ^{(\pm )}(\delta 𝚪,t)`$, namely

$`X(\delta 𝚪)_t^{(\pm )}{\displaystyle 𝑑\delta 𝚪X(\delta 𝚪)\rho ^{(\pm )}(\delta 𝚪,t)}`$ (35)
for any function $`X(\delta 𝚪)`$ of $`\delta 𝚪`$. In Ref.  the positive Lyapunov exponents were calculated by the time-averaged exponential rate of the diagonal elements of an orthogonal-transformed matrix of $`\mathrm{{\rm Y}}^{(+)}(t)`$. We get the negative Lyapunov exponents from the matrix $`\mathrm{{\rm Y}}^{()}(t)`$ by a similar procedure. The introduction of Lyapunov exponents by the spatial coordinate part only (or the momentum part only) of the tangent vector has been used previously by Refs. .

We introduce the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ defined by

$`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)\mathrm{{\rm Y}}^{(\pm )}(\pm t)e^{\pm \overline{\alpha }t}.`$ (36)
As shown in Appendix B, the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ satisfies the differential equation

$`{\displaystyle \frac{d^4\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^4}}{\displaystyle \frac{1}{2}}\left[\mathrm{\Omega }^2{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}+{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}\left(\mathrm{\Omega }^2\right)^T\right]`$ (37)
$`+{\displaystyle \frac{1}{16}}[\mathrm{\Omega }^4\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)2\mathrm{\Omega }^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)\left(\mathrm{\Omega }^2\right)^T`$ (38)
$`+\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)\left(\mathrm{\Omega }^4\right)^T]{\displaystyle \frac{2}{m^2}}\widehat{𝒟}\{{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}\}`$ (39)
$`=0`$ (40)
where we assumed the probability density $`\rho ^{(\pm )}(\delta 𝚪,t)`$ to be zero at the boundary of the tangent space. Here $`\mathrm{\Omega }`$ is the $`(2N)\times (2N)`$ matrix defined by

$`\mathrm{\Omega }\mathrm{\Phi }\mathrm{\Psi }.`$ (41)
and $`\widehat{𝒟}`$ is the linear operator to map any $`(2N)\times (2N)`$ matrix $`X(X_{jk})`$ to the $`(2N)\times (2N)`$ matrix $`\widehat{𝒟}\{X\}((\widehat{𝒟}\{X\})_{jk})`$ defined by

$`(\widehat{𝒟}\{X\})_{jk}{\displaystyle \underset{\mu =1}{\overset{N}{}}}{\displaystyle \underset{\nu =1}{\overset{N}{}}}D_{j\mu k\nu }X_{\mu \nu }.`$ (42)
It may be noted that Eq. (40) for the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ is invariant under the transformation $`\mathrm{\Omega }\mathrm{\Omega }`$.

We can choose the initial probability density $`\rho ^{(+)}(\delta 𝚪,0)`$ arbitrarily to calculate the positive Lyapunov exponents. On the other hand in order to derive the corresponding negative Lyapunov exponents we assume the initial probability density $`\rho ^{()}(\delta 𝚪,0)`$ for the time-reversed tangent vector to satisfy the condition $`d^k\stackrel{~}{\mathrm{{\rm Y}}}^{()}(t)/dt^k|_{t=0}=d^k\stackrel{~}{\mathrm{{\rm Y}}}^{(+)}(t)/dt^k|_{t=0}`$, $`k=0,1,2,3`$ at the initial time $`t=0`$. Under this assumption it follows from Eq. (40) that

$`\mathrm{{\rm Y}}^{()}(t)e^{\overline{\alpha }t}=\mathrm{{\rm Y}}^{(+)}(t)e^{\overline{\alpha }t},`$ (43)
because the quantities $`\stackrel{~}{\mathrm{{\rm Y}}}^{(+)}(t)`$ and $`\stackrel{~}{\mathrm{{\rm Y}}}^{()}(t)`$ defined by Eq. (36) satisfy the same differential equation (40) and have the same initial condition. Therefore the diagonal element $`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)`$ of any orthogonal-transformed matrix of $`\mathrm{{\rm Y}}^{(\pm )}(t)`$ must satisfy the relation

$`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{()}(t)=\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(+)}(t)e^{2\overline{\alpha }t}.`$ (44)
This equation connects the time-forward evolution and the time-reversed evolution in the amplitudes of the spatial part of the tangent vector in the thermostatted system.

The $`j`$-th positive (or zero) Lyapunov exponent $`\lambda _j^{(+)}`$ and its conjugate negative (or zero) exponent $`\lambda _j^{()}`$ are given by

$`\lambda _j^{(\pm )}=\underset{t=\pm \mathrm{}}{lim}{\displaystyle \frac{1}{2t}}\mathrm{ln}{\displaystyle \frac{\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)}{\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(0)}}.`$ (45)
The quantity $`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)`$ satisfies the condition $`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(+)}(0)=\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{()}(0)`$ at the initial time, and using Eqs. (44) and (45) we obtain

$`\lambda _j^{(+)}+\lambda _j^{()}=\overline{\alpha }.`$ (46)
This is the conjugate pairing rule of the Lyapunov spectrum for the iso-kinetic thermostatted system with a shear field. It is clear that it is attributed to the pairing rule of the Hamiltonian system in the case of $`\overline{\alpha }=0`$.

## V Conjugate pairing rule for a color Field

For an actual calculation of the Lyapunov spectrum for an iso-kinetic thermostatted system by the master equation approach we have to know the value of tensor $`D_{jkln}`$ and to solve the differential equation (40) for the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$. Let’s discuss these points briefly by using a case without a shear field such as a color field system , namely

$`\gamma =0.`$ (47)
For simplicity in this section we also assume that the tensor $`D_{jkln}`$ is expressed as the multiplication of the matrix elements of a symmetric $`(2N)\times (2N)`$ matrix $`W(W_{jk})`$:

$`D_{jkln}=W_{jk}W_{ln}.`$ (48)
This assumption has already been used in Ref. .

As shown in Appendix C, under the assumptions (47) and (48) the equation of the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ is simplified to

$`{\displaystyle \frac{d^3\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^3}}\overline{\alpha }^2{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}{\displaystyle \frac{2}{m^2}}W\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)W=0.`$ (49)
It is shown that Eq. (40) is given by taking the time-differential in both the sides of Eq. (49). Using the orthogonal matrix $`V`$ diagonalizing the matrix $`W`$, namely

$`(V^TWV)_{jk}=\omega _j\delta _{jk}`$ (50)
with a real eigenvalue $`\omega _j`$, the quantity $`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)`$ is expressed as the diagonal element of the matrix $`\stackrel{´}{\mathrm{{\rm Y}}}^{(\pm )}(t)(\stackrel{´}{\mathrm{{\rm Y}}}_{jk}^{(\pm )}(t))`$ defined by

$`\stackrel{´}{\mathrm{{\rm Y}}}^{(\pm )}(t)=V^T\mathrm{{\rm Y}}^{(\pm )}(t)V.`$ (51)
Using Eq. (51) we can solve the equation for the quantity $`\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)`$ derived from Eq. (49), and by using Eqs. (36), (45) and (51) the Lyapunov exponents are given by

$`\lambda _j^{(\pm )}={\displaystyle \frac{\overline{\alpha }}{2}}\pm {\displaystyle \frac{1}{2}}\left(\mathrm{\Lambda }_j+{\displaystyle \frac{\overline{\alpha }^2}{3\mathrm{\Lambda }_j}}\right).`$ (52)
where $`\mathrm{\Lambda }_j`$ is defined by

$`\mathrm{\Lambda }_j\left[\left({\displaystyle \frac{\omega _j}{m}}\right)^2+\sqrt{\left({\displaystyle \frac{\omega _j}{m}}\right)^4\left({\displaystyle \frac{\overline{\alpha }^2}{3}}\right)^3}\right]^{1/3}.`$ (53)
(See Appendix C about a derivation of Eq. (52).) It is clear that the Lyapunov exponents given by Eq. (52) satisfy the conjugate pairing rule (46).

Concerning the expression (52) for the Lyapunov exponent it is important to note that the tensor $`D_{jkln}`$ can depend on external force fields. This implies that the eigenvalue $`\omega _j`$ of the matrix $`W`$ can depend on the friction coefficient $`\overline{\alpha }`$. If we were to assume the quantity $`\omega _j`$ to be independent of the friction coefficient $`\overline{\alpha }`$, then we obtain the expression of the Lyapunov exponents as $`\lambda _j^{(\pm )}=\pm |\omega _j/(2m)|^{2/3}\overline{\alpha }/2+𝒪(\overline{\alpha }^2)`$ from Eq. (52) in the case of $`|\overline{\alpha }|\sqrt{3}|\omega _j/m|^{2/3}`$. However this is not consistent with the numerical results in a deterministic many-hard-disk system with a color field in which the negative Lyapunov exponents rather increase as the value of the friction coefficient increases . This consideration suggests that we should not neglect the external force field dependence of the correlation amplitude $`D_{jkln}`$ at least in the color field case. The dependence of the tensor $`D_{jkln}`$ on the shear rate and the external force fields should be a subject for a separated paper, although the conjugate pairing rule of the Lyapunov spectrum is correct regardless of their dependence as shown in this paper.

## VI Conclusion and Remarks

In this paper we have applied the master equation approach to Lyapunov spectra to non-equilibrium iso-kinetic thermostatted systems in order to discuss a conjugate pairing rule. We considered two-dimensional many-particle system with Gaussian white random interactions between the particles. In this system the positive Lyapunov exponents are calculated by a (forward) Fokker-Planck equation for the tangent vector dynamics. We proposed a method to calculate the negative Lyapunov exponents by a time-reversed master equation, especially the anti-Fokker-Planck equation where the diffusion term has the opposite sign to the forward Fokker-Planck equation. Using the Lyapunov exponents calculated by these two Fokker-Planck equations we showed the conjugate pairing rule of the Lyapunov spectrum for iso-kinetic thermostatted systems with a shear field given by the Sllod equation in the thermodynamic limit. We also gave an concrete form to connect the Lyapunov exponents with the time-correlation of the interaction matrix in a thermostatted system without a shear field.

We discussed the conjugate pairing rule based on the iso-kinetic thermostat in the thermodynamic limit. However it is known that the iso-kinetic thermostat is formally equivalent to other thermostats such as the iso-energetic thermostat in the thermodynamic limit . In this sense our result should be correct in systems with such other thermostats, more explicitly as far as the friction coefficient can be regarded as a constant even in a finite number of particle systems.

In order to get the ani-Fokker-Planck equation we used the fact that the shear rate $`\gamma `$ changes its sign in the time reversed motion. However this time-reversed change of the sign of the shear rate to get the anti-Fokker-Planck equation may not be essential to obtain the negative Lyapunov exponents, if the Lyapunov exponents are even functions of the shear rate. We have not proved the invariance of the Lyapunov exponents under the transformation $`\gamma \gamma `$ in this paper, but the invariance of Eq. (40) under the transformation $`\mathrm{\Omega }\mathrm{\Omega }`$ implies that the Lyapunov exponents are invariant under this transformation.

We can show that all the Lyapunov exponents $`\lambda _j^{(+)}`$ ($`\lambda _j^{()}`$) are non-negative (non-positive) in the case of $`\gamma =0`$ (See Appendix C.). This implies that in this case the number of the positive Lyapunov exponents should be equal to the number of the negative Lyapunov exponents. However we have not proved that it is also correct in the presence of a shear field: $`\gamma 0`$. Concerning this point we should notice that a numerical calculation of the Lyapunov spectrum for the Sllod equation (1) and (2) showed that in the case of a high shear rate the number of the positive Lyapunov exponents can be less than the number of the negative Lyapunov exponents . Therefore it should be interesting to check whether the master equation approach to the Lyapunov spectrum can describe such a situation or not.

It should be noted that the discussion of this paper does not conclude that the conjugate pairing rule of the Lyapunov spectrum for the thermostatted system with a shear field must be satisfied rigorously in deterministic chaotic systems. To show the conjugate pairing rule in this paper we assumed the Gaussian white randomness (14) and (18) for the particle interactions, and there is no guarantee that we can justify the conjugate pairing rule by the master equation approach under a more general random interaction of particles, especially under the randomness with a long time correlation which leads to a more general master equation for the tangent vector rather than a simple Fokker-Planck equation. A generalization of the conjugate pairing rule by the master equation approach to a more general random interaction remains as one of the important future problems.

## Acknowledgements

We are grateful for financial support for this work from the Australian Research Council.

## A Master Equation for the Tangent Vector

In this appendix we derive the Fokker-Planck equation (23) for the tangent vector space. Using the Kramers-Moyal expansion the dynamics of the probability density $`\rho ^{(+)}(\delta 𝚪,t)`$ is given by

$`{\displaystyle \frac{\rho ^{(+)}(\delta 𝚪,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{j_1=1}{\overset{2N}{}}}{\displaystyle \underset{j_2=1}{\overset{2N}{}}}\mathrm{}{\displaystyle \underset{j_n=1}{\overset{2N}{}}}(1)^n`$ (A.2)
$`\times {\displaystyle \frac{^n\mathrm{\Xi }_{j_1j_2\mathrm{}j_n}^{(n)}(\delta 𝚪,t)\rho ^{(+)}(\delta 𝚪,t)}{\delta \Gamma _{j_1}\delta \Gamma _{j_2}\mathrm{}\delta \Gamma _{j_n}}}`$
where $`\mathrm{\Xi }_{j_1j_2\mathrm{}j_n}^{(n)}(\delta 𝚪,t)`$ is defined by

$`\mathrm{\Xi }_{j_1j_2\mathrm{}j_n}^{(n)}(\delta 𝚪,t)`$ (A.3)
$`{\displaystyle \frac{1}{n!}}\underset{s0}{lim}{\displaystyle \frac{1}{s}}[\delta \Gamma _{j_1}(t+s)\delta \Gamma _{j_1}(t)]`$ (A.4)
$`\times [\delta \Gamma _{j_2}(t+s)\delta \Gamma _{j_2}(t)]`$ (A.5)
$`\mathrm{}\times [\delta \Gamma _{j_n}(t+s)\delta \Gamma _{j_n}(t)]|_{\delta 𝚪(t)=\delta 𝚪}`$ (A.6)
and $`\delta \Gamma _j(t)`$ is the $`j`$-th component of the tangent vector $`\delta 𝚪(t)`$.

Using Eq. (7) we obtain

$`\delta 𝚪(t+s)\delta 𝚪(t)`$ (A.7)
$`=\{\stackrel{}{T}\mathrm{exp}\left[{\displaystyle _t^{t+s}}d\tau (\tau )\right]1\}\delta 𝚪(t)`$ (A.8)
$`={\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle _t^{t+s}}𝑑\tau _n{\displaystyle _t^{\tau _n}}𝑑\tau _{n1}\mathrm{}{\displaystyle _t^{\tau _2}}𝑑\tau _1`$ (A.9)
$`\times (\tau _n)(\tau _{n1})\mathrm{}(\tau _1)\delta 𝚪(t).`$ (A.10)
It follows from Eqs. (10), (14), (18), (A.6) and (A.10) that

$`𝚵^{(1)}(\delta 𝚪,t)`$ (A.11)
$`(\mathrm{\Xi }_1^{(1)}(\delta 𝚪,t),\mathrm{\Xi }_2^{(1)}(\delta 𝚪,t),\mathrm{},\mathrm{\Xi }_{2N}^{(1)}(\delta 𝚪,t))^T`$ (A.12)
$`=\underset{s0}{lim}{\displaystyle \frac{1}{s}}\left[\delta 𝚪(t+s)\delta 𝚪(t)\right]|_{\delta 𝚪(t)=\delta 𝚪}`$ (A.13)
$`=\underset{s0}{lim}{\displaystyle \frac{1}{s}}{\displaystyle _t^{t+s}}𝑑\tau (\tau )\delta 𝚪`$ (A.14)
$`=\left(\begin{array}{c}\mathrm{\Phi }\delta 𝐪+\delta 𝐩/m\\ \mathrm{\Psi }\delta 𝐩/m\end{array}\right)`$ (A.17)
and

$`\mathrm{\Xi }^{(2)}(\delta 𝚪,t)(\mathrm{\Xi }_{jk}^{(2)}(\delta 𝚪,t))`$ (A.18)
$`=\underset{s0}{lim}{\displaystyle \frac{1}{2s}}[\delta 𝚪(t+s)\delta 𝚪(t)]`$ (A.19)
$`\times [\delta 𝚪(t+s)\delta 𝚪(t)]^T|_{\delta 𝚪(t)=\delta 𝚪}`$ (A.20)
$`=\underset{s0}{lim}{\displaystyle \frac{1}{2s}}{\displaystyle _t^{t+s}}𝑑\kappa {\displaystyle _t^{t+s}}𝑑\tau `$ (A.21)
$`\times (\kappa )\delta 𝚪\delta 𝚪^T(\tau )^T`$ (A.22)
$`=\left(\begin{array}{cc}0_{2N}& 0_{2N}\\ 0_{2N}& \eta (\delta 𝐪)\end{array}\right)`$ (A.25)
where $`\eta (\delta 𝐪)(\eta _{jk}(\delta 𝐪))`$ is defined by

$`\eta _{jk}(\delta 𝐪){\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu =1}{\overset{N}{}}}{\displaystyle \underset{\nu =1}{\overset{N}{}}}D_{j\mu k\nu }\delta q_\mu \delta q_\nu .`$ (A.26)
Here the only non-zero contributions come from the $`n=1`$ term of Eq. (A.10). For general $`n`$, the number of delta functions from Eq. (18) must be only one less than the number of time integrals, to give a non-zero contribution. It is straightforward to show that this never happens for $`n>1`$. Concerning the terms including $`\mathrm{\Xi }_{j_1j_2\mathrm{}j_n}^{(n)}(\delta 𝚪,t)`$, $`n=3,4,\mathrm{}`$ in the right-hand side of Eq. (A.2) we obtain

$`\mathrm{\Xi }_{j_1j_2\mathrm{}j_n}^{(n)}(\delta 𝚪,t)=0,n=3,4,\mathrm{},`$ (A.27)
because of the Gaussian white properties (14) and (18) of the random matrix $`R(t)`$. Using Eqs. (A.2), (A.17), (A.25) and (A.27) we obtain the Fokker-Planck equation (23).

## B Equation for the Matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}`$

In this appendix we give details of the derivation of Eq. (40) from Eqs. (23), (33) and (34). We start this derivation by introducing the $`N\times N`$ matrices $`^{(\pm )}(t)`$ and $`𝒢^{(\pm )}(t)`$ defined by

$`^{(\pm )}(t)\delta 𝐪\delta 𝐩^T_t^{(\pm )},`$ (B.1)
$`𝒢^{(\pm )}(t)\delta 𝐩\delta 𝐩^T_t^{(\pm )}.`$ (B.2)
Equations (23) and (33) lead to the equations

$`{\displaystyle \frac{d\mathrm{{\rm Y}}^{(\pm )}(t)}{dt}}`$ $`=`$ $`\mathrm{\Phi }\mathrm{{\rm Y}}^{(\pm )}(t)+\mathrm{{\rm Y}}^{(\pm )}(t)\mathrm{\Phi }^T`$ (B.4)
$`+{\displaystyle \frac{1}{m}}[^{(\pm )}(t)+^{(\pm )}(t){}_{}{}^{T}],`$
$`{\displaystyle \frac{d^{(\pm )}(t)}{dt}}`$ $`=`$ $`\mathrm{\Phi }^{(\pm )}(t)+^{(\pm )}(t)\mathrm{\Psi }^T+{\displaystyle \frac{1}{m}}𝒢^{(\pm )}(t),`$ (B.5)
$`{\displaystyle \frac{d𝒢^{(\pm )}(t)}{dt}}`$ $`=`$ $`\mathrm{\Psi }𝒢^{(\pm )}(t)+𝒢^{(\pm )}(t)\mathrm{\Psi }^T\pm \widehat{𝒟}\left\{\mathrm{{\rm Y}}^{(\pm )}(t)\right\}`$ (B.7)
for the matrices $`\mathrm{{\rm Y}}^{(\pm )}(t)`$, $`^{(\pm )}(t)`$ and $`𝒢^{(\pm )}(t)`$ with the operator $`\widehat{𝒟}`$ defined by Eq. (42). Here, to derive Eq. (LABEL:ConneCalG) we used the relation

$`D_{lnjk}=D_{jkln}`$ (B.9)
which is derived from the definition (18) of the tensor $`D_{jkln}`$ and the symmetry property of the matrix $`R(t)`$.

Eqs. (B.4), (LABEL:ConneCalF) and (LABEL:ConneCalG) are equivalent to

$`{\displaystyle \frac{d\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}=\stackrel{˘}{}^{(\pm )}(t)P(t)^T+P(t)\stackrel{˘}{}^{(\pm )}(t){}_{}{}^{T},`$ (B.10)
$`{\displaystyle \frac{d\stackrel{˘}{}^{(\pm )}(t)}{dt}}=P(t)\stackrel{˘}{𝒢}^{(\pm )}(t),`$ (B.11)
$`{\displaystyle \frac{d\stackrel{˘}{𝒢}^{(\pm )}(t)}{dt}}=\pm \widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}`$ (B.12)
where $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$, $`\stackrel{˘}{}^{(\pm )}(t)`$ and $`\stackrel{˘}{𝒢}^{(\pm )}(t)`$ are defined by

$`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ $``$ $`e^{\mathrm{\Phi }t}\mathrm{{\rm Y}}^{(\pm )}(t)e^{\mathrm{\Phi }^Tt},`$ (B.13)
$`\stackrel{˘}{}^{(\pm )}(t)`$ $``$ $`e^{\mathrm{\Phi }t}^{(\pm )}(t)e^{\mathrm{\Psi }^Tt},`$ (B.14)
$`\stackrel{˘}{𝒢}^{(\pm )}(t)`$ $``$ $`e^{\mathrm{\Psi }t}𝒢^{(\pm )}(t)e^{\mathrm{\Psi }^Tt},`$ (B.15)
and $`P(t)`$ is the $`(2N)\times (2N)`$ matrix defined by

$`P(t){\displaystyle \frac{1}{m}}e^{\mathrm{\Omega }t}`$ (B.16)
with the matrix $`\mathrm{\Omega }`$ defined by Eq. (41), and $`\widehat{\stackrel{~}{D}}_t\{\mathrm{}\}`$ is defined by

$`\widehat{\stackrel{~}{𝒟}}_t\{X\}e^{\mathrm{\Psi }t}\widehat{𝒟}\left\{e^{\mathrm{\Phi }t}Xe^{\mathrm{\Phi }^Tt}\right\}e^{\mathrm{\Psi }^Tt}`$ (B.17)
for any $`(2N)\times (2N)`$ matrix $`X`$. Here we used the relation

$`\mathrm{\Phi }\mathrm{\Psi }=\mathrm{\Psi }\mathrm{\Phi },`$ (B.18)
so that we have the equation $`\mathrm{exp}\{\mathrm{\Phi }t\}\mathrm{exp}\{\mathrm{\Psi }t\}=\mathrm{exp}\{(\mathrm{\Phi }\mathrm{\Psi })t\}`$.

Noting that the matrix $`\stackrel{˘}{𝒢}^{(\pm )}(t)`$ is symmetric and the inverse matrix of the matrix $`P(t)`$ is given by $`P(t)^1=m\mathrm{exp}\{\mathrm{\Omega }t\}`$, we obtain

$`2\stackrel{˘}{𝒢}^{(\pm )}(t)`$ $`=`$ $`P(t)^1{\displaystyle \frac{d\stackrel{˘}{}^{(\pm )}(t)}{dt}}+{\displaystyle \frac{d\stackrel{˘}{}^{(\pm )}(t)^T}{dt}}\left[P(t)^1\right]^T`$ (B.19)
$`=`$ $`P(t)^1{\displaystyle \frac{d^2\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}\left[P(t)^1\right]^T`$ (B.22)
$`+P(t)^1\stackrel{˘}{}^{(\pm )}(t)\mathrm{\Omega }^T`$
$`+\mathrm{\Omega }\stackrel{˘}{}^{(\pm )}(t){}_{}{}^{T}\left[P(t)^1\right]_{}^{T}`$
by using Eqs. (B.10), (B.11) and (B.16). Besides, using Eqs. (B.10), (B.11) and (B.16) we obtain

$`{\displaystyle \frac{d}{dt}}\{P(t)^1\stackrel{˘}{}^{(\pm )}(t)\mathrm{\Omega }^T+\mathrm{\Omega }\stackrel{˘}{}^{(\pm )}(t){}_{}{}^{T}\left[P(t)^1\right]_{}^{T}\}`$ (B.23)
$`=\mathrm{\Omega }P(t)^1{\displaystyle \frac{d\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}\left[P(t)^1\right]^T\mathrm{\Omega }^T`$ (B.24)
$`+\stackrel{˘}{𝒢}^{(\pm )}(t)\mathrm{\Omega }^T+\mathrm{\Omega }\stackrel{˘}{𝒢}^{(\pm )}(t)`$ (B.25)
where we again used the relation $`\stackrel{˘}{𝒢}^{(\pm )}(t)^T=\stackrel{˘}{𝒢}^{(\pm )}(t)`$. It follows from Eqs. (B.10), (B.11) and (B.12) that

$`\pm 2\widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}`$ (B.26)
$`={\displaystyle \frac{d}{dt}}P(t)^1{\displaystyle \frac{d^2\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}\left[P(t)^1\right]^T`$ (B.27)
$`+\mathrm{\Omega }P(t)^1{\displaystyle \frac{d\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}\left[P(t)^1\right]^T\mathrm{\Omega }^T`$ (B.28)
$`+\stackrel{˘}{𝒢}^{(\pm )}(t)\mathrm{\Omega }^T+\mathrm{\Omega }\stackrel{˘}{𝒢}^{(\pm )}(t).`$ (B.29)
Taking the time-differential of both the sides of Eq. (B.29), and using Eq. (B.12) we obtain

$`{\displaystyle \frac{d^2}{dt^2}}P(t)^1{\displaystyle \frac{d^2\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}\left[P(t)^1\right]^T`$ (B.30)
$`+\mathrm{\Omega }{\displaystyle \frac{d}{dt}}P(t)^1{\displaystyle \frac{d\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}\left[P(t)^1\right]^T\mathrm{\Omega }^T`$ (B.31)
$`2{\displaystyle \frac{d}{dt}}\widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}`$ (B.32)
$`\pm \widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}\mathrm{\Omega }^T\pm \mathrm{\Omega }\widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}`$ (B.33)
$`=0.`$ (B.34)
This is the equation for the quantity $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ only.

Now we derive the equation for $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ defined by Eq. (36) from Eq. (B.34) for $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ defined by Eq. (B.13). We note

$`\mathrm{\Omega }`$ $`=`$ $`2\mathrm{\Phi }+\overline{\alpha }I_{2N}`$ (B.35)
$`=`$ $`2\mathrm{\Psi }\overline{\alpha }I_{2N}`$ (B.36)
which is derived from Eqs. (11), (12) and (41). Using Eqs. (36), (B.13) and (B.35) the matrices $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ are connected with the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ by

$`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)=e^{\mathrm{\Omega }t/2}\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)e^{\mathrm{\Omega }^Tt/2}.`$ (B.37)
We introduce the multiplication $`XY`$ of $`X`$ and $`Y`$ which is defined by

$`XY{\displaystyle \frac{1}{2}}\left[XY+(XY)^T\right]`$ (B.38)
for any square matrices $`X`$ and $`Y`$ of the same size. This multiplication is used in the relation

$`{\displaystyle \frac{d}{dt}}e^{\pm \mathrm{\Omega }t/2}Z(t)e^{\pm \mathrm{\Omega }^Tt/2}`$ (B.39)
$`=e^{\pm \mathrm{\Omega }t/2}\left[{\displaystyle \frac{dZ(t)}{dt}}\pm \mathrm{\Omega }Z(t)\right]e^{\pm \mathrm{\Omega }^Tt/2}`$ (B.40)
satisfied by any $`(2N)\times (2N)`$ symmetric matrix $`Z(t)`$ as a function of $`t`$. Noting Eqs. (B.38) and (B.40), Eq. (B.37) leads to

$`P(t)^1{\displaystyle \frac{d\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt}}\left[P(t)^1\right]^T`$ (B.41)
$`=m^2e^{\mathrm{\Omega }t/2}\left[{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt}}\mathrm{\Omega }\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\right]e^{\mathrm{\Omega }^Tt/2},`$ (B.42)
$`P(t)^1{\displaystyle \frac{d^2\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)}{dt^2}}\left[P(t)^1\right]^T`$ (B.43)
$`=m^2e^{\mathrm{\Omega }t/2}\{{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^2}}2\mathrm{\Omega }{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt}}`$ (B.44)
$`+\mathrm{\Omega }[\mathrm{\Omega }\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)]\}e^{\mathrm{\Omega }^Tt/2}`$ (B.45)
where we used the relation $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)^T=\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$. Moreover $`\widehat{\stackrel{~}{𝒟}}_t`$ operated on the matrix $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ is connected with $`\widehat{𝒟}`$ operated on the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)`$ as

$`\widehat{\stackrel{~}{𝒟}}_t\left\{\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)\right\}=e^{\mathrm{\Omega }t/2}\widehat{𝒟}\left\{\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\right\}e^{\mathrm{\Omega }^Tt/2}`$ (B.46)
where we used Eqs. (36), (B.13) (B.17) and (B.36). Inserting Eqs. (B.42), (B.45) and (B.46) into Eq. (B.34) and using Eq. (B.40) we obtain

$`{\displaystyle \frac{d^4\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^4}}2\mathrm{\Omega }\left[\mathrm{\Omega }{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^2}}\right]`$ (B.47)
$`+\mathrm{\Omega }{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^2}}\mathrm{\Omega }^T`$ (B.48)
$`+\mathrm{\Omega }\left[\mathrm{\Omega }\left[\mathrm{\Omega }\left[\mathrm{\Omega }\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\right]\right]\right]`$ (B.49)
$`\mathrm{\Omega }\left[\mathrm{\Omega }\left[\mathrm{\Omega }\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\right]\right]\mathrm{\Omega }^T`$ (B.50)
$`{\displaystyle \frac{2}{m^2}}\widehat{𝒟}\left\{{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt}}\right\}=\mathrm{\hspace{0.33em}\hspace{0.33em}0}.`$ (B.51)
Eq. (B.51) is equivalent to

$`{\displaystyle \frac{d^4\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^4}}`$ (B.52)
$`{\displaystyle \frac{1}{2}}\left[\mathrm{\Omega }^2{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^2}}+{\displaystyle \frac{d^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^2}}\left(\mathrm{\Omega }^2\right)^T\right]`$ (B.53)
$`+{\displaystyle \frac{1}{16}}[\mathrm{\Omega }^4\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)2\mathrm{\Omega }^2\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\left(\mathrm{\Omega }^2\right)^T`$ (B.54)
$`+\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\left(\mathrm{\Omega }^4\right)^T]{\displaystyle \frac{2}{m^2}}\widehat{𝒟}\{{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt}}\}`$ (B.55)
$`=0.`$ (B.56)
By exchanging $`t`$ with $`\pm t`$ in Eq. (B.56) we obtain Eq. (40).

## C Lyapunov Exponents in the color Field Case

In this appendix we consider the case with no shear field using the condition (47), and derive Eq. (49) under the assumption (48). We also give a derivation of Eq. (52) briefly.

Under the condition (47), the matrix $`\mathrm{\Omega }`$ defined by Eq. (41) is simply an identical matrix multiplied by a constant and is given by

$`\mathrm{\Omega }=\overline{\alpha }I_{2N},`$ (C.1)
and the matrix $`P(t)`$ defined by Eq. (B.16) and the operator $`\widehat{\stackrel{~}{𝒟}}_t\{\mathrm{}\}`$ defined by Eq. (B.17) are given by

$`P(t)={\displaystyle \frac{1}{m}}e^{\overline{\alpha }t}I_{2N},`$ (C.2)
$`\widehat{\stackrel{~}{𝒟}}_t\left\{X\right\}=e^{2\overline{\alpha }t}\widehat{𝒟}\left\{X\right\}`$ (C.3)
for any $`(2N)\times (2N)`$ matrix $`X`$. Noting Eqs. (C.2) and (C.3) and the relation $`\stackrel{˘}{\mathrm{{\rm Y}}}^{(\pm )}(t)=\mathrm{{\rm Y}}^{(\pm )}(t)`$, Eqs. (B.10), (B.11) and (B.12) are simply attributed into

$`{\displaystyle \frac{d\mathrm{{\rm Y}}^{(\pm )}(t)}{dt}}={\displaystyle \frac{1}{m}}e^{\overline{\alpha }t}[\stackrel{˘}{}^{(\pm )}(t)+\stackrel{˘}{}^{(\pm )}(t){}_{}{}^{T}],`$ (C.4)
$`{\displaystyle \frac{d\stackrel{˘}{}^{(\pm )}(t)}{dt}}={\displaystyle \frac{1}{m}}e^{\overline{\alpha }t}\stackrel{˘}{𝒢}^{(\pm )}(t),`$ (C.5)
$`{\displaystyle \frac{d\stackrel{˘}{𝒢}^{(\pm )}(t)}{dt}}=\pm e^{2\overline{\alpha }t}\widehat{𝒟}\left\{\mathrm{{\rm Y}}^{(\pm )}(t)\right\}.`$ (C.6)
Noting that the matrix $`\stackrel{˘}{𝒢}^{(\pm )}(t)`$ is symmetric, Eqs. (C.4), (C.5) and (C.6) lead to the differential equation

$`{\displaystyle \frac{d}{dt}}e^{\overline{\alpha }t}{\displaystyle \frac{d}{dt}}e^{\overline{\alpha }t}{\displaystyle \frac{d}{dt}}\mathrm{{\rm Y}}^{(\pm )}(t)=\pm {\displaystyle \frac{2}{m^2}}e^{2\overline{\alpha }t}\widehat{𝒟}\left\{\mathrm{{\rm Y}}^{(\pm )}(t)\right\}`$ (C.7)
for the matrix $`\mathrm{{\rm Y}}^{(\pm )}(t)`$ only.

Now we consider the derivation of the equation for the matrix $`\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(t)`$ (defined by Eq. (36)) from Eq. (C.7) for the matrix $`\mathrm{{\rm Y}}^{(\pm )}(t)`$. It follows from Eqs. (36) and (C.7) that

$`{\displaystyle \frac{d^3\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt^3}}\overline{\alpha }^2{\displaystyle \frac{d\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)}{dt}}{\displaystyle \frac{2}{m^2}}\widehat{𝒟}\left\{\stackrel{~}{\mathrm{{\rm Y}}}^{(\pm )}(\pm t)\right\}=0.`$ (C.8)
(C.9)
By exchanging $`t`$ with $`\pm t`$ in Eq. (C.9) and using Eqs. (42) and (48) we obtain Eq. (49).

Using Eqs. (49) and (51) the real function $`\xi _j^{(\pm )}(t)`$ of $`t`$ defined by

$`\xi _j^{(\pm )}(t)(V^T\stackrel{~}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(t)V)_{jj}=\stackrel{´}{\mathrm{{\rm Y}}}_{jj}^{(\pm )}(\pm t)e^{\pm \overline{\alpha }t}`$ (C.10)
satisfies the equation

$`{\displaystyle \frac{d^3\xi _j^{(\pm )}(t)}{dt^3}}\overline{\alpha }^2{\displaystyle \frac{d\xi _j^{(\pm )}(t)}{dt}}{\displaystyle \frac{2\omega _j^2}{m^2}}\xi _j^{(\pm )}(t)=0.`$ (C.11)
The real solution of the linear differential equation (C.11) is expressed as

$`\xi _j^{(\pm )}(t)={\displaystyle \underset{k=1}{\overset{3}{}}}\text{Re}\left\{\upsilon _j^{(k)}e^{\zeta _j^{(k)}t}\right\}`$ (C.12)
where $`\upsilon _j^{(k)}`$, $`j=1,2,3`$ are constants determined by the initial condition, and $`\zeta _j^{(k)}`$, $`j=1,2,3`$ are the three solutions of the equation

$`\zeta ^3\overline{\alpha }^2\zeta {\displaystyle \frac{2\omega _j^2}{m^2}}=0`$ (C.13)
for $`\zeta `$. Here $`\text{Re}\{X\}`$ means to take the real part of any imaginary number $`X`$. We sort the quantities $`\zeta _j`$, $`j=1,2,3`$ as $`\text{Re}\{\zeta _j^{(1)}\}\text{Re}\{\zeta _j^{(2)}\}\text{Re}\{\zeta _j^{(3)}\}`$, so that using Eqs. (45), (51) and (C.10) the Lyapunov exponent $`\lambda _j^{(\pm )}`$ is expressed as

$`\lambda _j^{(\pm )}`$ $`=`$ $`\pm \underset{t=+\mathrm{}}{lim}{\displaystyle \frac{1}{2t}}\mathrm{ln}{\displaystyle \frac{\xi _j^{(\pm )}(t)e^{\overline{\alpha }t}}{\xi _j^{(\pm )}(0)}}`$ (C.14)
$`=`$ $`{\displaystyle \frac{\overline{\alpha }}{2}}\pm {\displaystyle \frac{1}{2}}\text{Re}\{\zeta _j^{(1)}\}.`$ (C.15)
It follows from Eq. (C.13) that $`\zeta =\zeta _j^{(1)}`$ is a real solution of Eq. (C.13) and satisfies the conditions $`\zeta _j^{(1)}|\overline{\alpha }|`$ $`lim_{\omega _j0}\zeta _j^{(1)}=|\overline{\alpha }|`$, noting that the point $`\zeta =\zeta _j^{(1)}`$ is the maximum intersecting point of the graphs $`y=\zeta ^3\overline{\alpha }^2\zeta `$ and $`y=2\omega _j^2/m^2`$ in the $`\zeta `$-$`y`$ plain. This means that the Lyapunov exponents $`\lambda _j^{(+)}`$ ($`\lambda _j^{()}`$) must be non-negative (non-positive). More concretely the quantity $`\zeta _j^{(1)}`$ is given by

$`\zeta _j^{(1)}=\mathrm{\Lambda }_j+{\displaystyle \frac{\overline{\alpha }^2}{3\mathrm{\Lambda }_j}}`$ (C.16)
with the quantity $`\mathrm{\Lambda }_j`$ defined by Eq. (53). We can check that in the case of $`|\overline{\alpha }|\sqrt{3}|\omega _j/m|^{2/3}`$ the quantities $`\mathrm{\Lambda }_j`$ and $`\zeta _j^{(1)}`$ are both real numbers, and in the case of $`|\overline{\alpha }|>\sqrt{3}|\omega _j/m|^{2/3}`$ the quantity $`\mathrm{\Lambda }_j`$ can be an imaginary number but the quantity $`\zeta _j^{(1)}`$ given by Eq. (C.16) is still a real number and satisfies the condition $`lim_{\omega _j0}\zeta _j^{(1)}=|\overline{\alpha }|`$. Using Eqs. (C.15) and (C.16) we obtain Eq. (52).