text:
# Perspectives on Finding the Neutrino Nature

\[

## Abstract

The possibility of determining the neutrino nature is considered in view of the most recent experimental observations. The analysis combines schemes with three and four neutrinos.The data on oscillations is put together with that from the search of neutrinoless double beta decay and results on tritium beta decay. All solar neutrino oscillation solutions are taken into account. The sensitivity of the problem on future experimental bounds from GENIUS is studied. Dirac neutrinos are shown to be unavoidable already at present in some schemes and the constraints will quickly become more stringent with future data. The consequences of including bounds from Cosmology on the neutrino content of Hot Dark Matter are commented.

\]

At present, particle physics seems to have a long awaited opportunity to discover the path towards an extension of the Standard Model. This opportunity is given by the recent confirmation of neutrino oscillations . With the currently accumulated experimental data, we have many hints on the size of neutrino masses and mixing angles. This information imposes serious constraints on the structure of the model. Unfortunately there remains a major unknown, which is the neutrino nature. As a neutral particle it can be described either by means of the Dirac equation or by means of the Majorana equation . This problem is of the utmost importance for model builders.

It has been proved many times, that with the current precision, the only experiment able to show which of the two descriptions is valid, is the search for the neutrinoless double beta decay $`(\beta \beta )_{0\nu }`$. Should this process be observed, we would be left with no other choice but to acknowledge neutrinos to be Majorana particles. As this has not been the case, we have at hand only a lower bound on the lifetime of the $`(\beta \beta )_{0\nu }`$ decaying nuclei. This alone is not enough. However when supplied with the results of neutrino oscillation experiments and the tritium beta decay mass bounds, interesting, but this time neutrino mass scheme dependent conclusions can be drawn. The situation will undoubtedly improve, when the future GENIUS  experiment will furnish its first results.

In this paper, we combine the neutrino oscillation data, the tritium beta decay bounds and the present and future GENIUS bounds on the $`(\beta \beta )_{0\nu }`$ lifetime (if their results are not positive, since in the opposite case the situation is clear, as stated above) in order to derive conditions which would define the neutrino as a Dirac particle. To this end we take into account the possibility of having three or four light neutrinos. Therefore our conclusions will be valid in either case, whether the LSND  results would eventually be ruled out or confirmed.

The novelty in respect to previous analyses  is the consideration of the case of three and four neutrinos together, supplemented with the most recent bounds on the mixing matrix elements and the variation of the minimal neutrino mass $`(m_\nu )_{min}`$ in its full range. Furthermore we do not assume that neutrinos are Majorana particles like in so many other works .

The oscillation data is described in terms of differences of neutrino mass squares $`\delta m_{ij}^2m_i^2m_j^2`$. Define a scheme of neutrino masses as an ordering of these $`\delta m^2`$’s. For three neutrinos there are two such schemes, whereas for four neutrinos there are six. In the latter case, four schemes can be rejected by means of a combined analysis . All four possibilities are shown in Fig. 2 and Fig 3. It should be clear that with the strong ordering of $`\delta m^2`$’s ($`\delta m_{sol}^2\delta m_{atm}^2\delta m_{LSND}^2`$), it is only the largest that actually matters. Therefore for three neutrinos we need only to know $`\delta m_{atm}^23\times 10^3eV^2`$ , whereas for four $`\delta m_{LSND}^21eV^2`$ . The second important information is that the spectrum of the masses of light neutrinos is bounded from above by $`m_\nu <2.5eV`$ . With the additional information on the value of one mass, we recover the complete mass spectrum of neutrinos.

The basic quantity to consider is the effective neutrino mass, measured in $`(\beta \beta )_{0\nu }`$ decay

$$|m_\nu |\left|\underset{i=light}{}U_{ei}^2m_i\right|,$$
(1)
with $`U`$ the mixing matrix which we assume unitary. Current and future bounds on $`|m_\nu |`$ are given in Table I. At present, nothing can be said about the phases of $`U`$, therefore we allow them to have the least favorable values. Our analysis is based on the following statement: if the minimum of $`|m_\nu |`$ with respect to $`U_{ei}`$ within their allowed range, $`(|m_\nu |)_{min}`$, exceeds the experimental bound, then neutrinos are Dirac particles .

Suppose that the solar neutrino oscillations take place mostly between the first and the second neutrino states, then

$$(|m_\nu |)_{min}=\mathrm{min}\left|l^em_as^em_b\right|,$$
(2)
with

$$l^e[l_{}^e,l_+^e],s^e[s_{}^e,s_+^e],$$
(3)
and

$$l_+^e=(1s_+^e),l_{}^e=(1s_+^e)\sqrt{1\mathrm{sin}^22\theta _{sun}},$$
(4)
where in the case of three neutrinos

$$s_+^e=s_{}^e=|U_{e3}|^2,$$
(5)
and in the case of four

$$s_+^e=|U_{e3}|^2+|U_{e4}|^2,s_{}^e=\left||U_{e3}|^2|U_{e4}|^2\right|.$$
(6)
The mass parameters in Eq. 2 can be chosen in two different ways, depending on the ordering of masses of the neutrino states.

Schemes $`A_3`$ and $`A_4`$. The behavior of the minimum of the effective mass is the following

$$(|m_\nu |)_{min}=\{\begin{array}{c}s_{}^e\sqrt{(m_\nu )_{min}^2+\delta m^2}l_+^e(m_\nu )_{min},\\ 0,\\ l_{}^e(m_\nu )_{min}s_+^e\sqrt{(m_\nu )_{min}^2+\delta m^2},\end{array}$$
(7)
with $`\delta m^2\delta m_{atm}^2`$ for three neutrinos and $`\delta m^2\delta m_{LSND}^2`$ for four. The third range occurs only if $`l_{}^e>s_+^e`$. This can be rewritten as

$$\mathrm{sin}^22\theta _{sun}<\frac{12s_+^e}{(1s_+^e)^2},$$
(8)
Obviously these schemes always allow for a Majorana mass. If however the condition Eq. 8 is satisfied, then it is possible to derive a lower bound on the mass spectrum which excludes this possibility. This bound will be larger than the value of $`(m_\nu )_{min}`$ which marks the beginning of the third range in Eq. 7

$$(m_\nu )_{min}=\left(\frac{\delta m^2(s_+^e)^2}{(l_{}^e)^2(s_+^e)^2}\right)^{1/2}.$$
(9)
This number will be small in the case of three neutrinos if $`\mathrm{sin}^22\theta _{sun}`$ is not too close to one, due to the small $`\delta m_{atm}^2`$, but can quickly grow large in the case of four.

Schemes $`B_3`$ and $`B_4`$. At the condition Eq. 8 we have

$$(|m_\nu |)_{min}=l_{}^e\sqrt{(m_\nu )_{min}^2+\delta m^2}s_+^e(m_\nu )_{min}.$$
(10)
Obviously, if the masses of the neutrinos are distributed according to these schemes the neutrino nature can be determined irrespective of the lowest mass, depending solely on the experimental bound on $`|m_\nu |`$

Dependence on the mixing matrix. The condition Eq. 8 plays a crucial role in the problem. The availability of lower bounds on the neutrino mass is tightly connected to it. We therefore have to investigate the size of both sides of the inequality.

The quantity $`s_+^e`$ can be fitted to various experimental data. It turns out that the solar neutrino flux observation is the least constraining since one has a bound of $`s_+^e<.67.73`$ . Remark however that the difference between $`\chi ^2`$ for $`s_+^e=0`$ and the above result is almost irrelevant. Therefore if other experiments constrain the considered quantity to be close to zero, then solar neutrino oscillation should be viewed as a two state phenomenon and 3 neutrino fits to solar data alone seem to be of little value. This will be important subsequently. The currently strongest constraint comes from a combined analysis of atmospheric and reactor data. Namely one has $`s_+^e<.01.02`$ . However a recent fit to the new (830-920 days) data of Super-Kamionkande  gives a best value of $`s_+^e.03`$. There has also been some speculation about strengthening this bound by adding some theoretical assumptions on the shape of the mass matrix. It seems that $`s_+^e0`$ is the most preferred value . In Fig. 1, the dependence of the upper bound on $`\mathrm{sin}^22\theta _{sun}`$ is shown as function of $`s_+^e`$. In the range $`s_+^e[0,.04]`$ this bound is between $`.998`$ and $`1`$.

Let us now turn to the experimental bounds on $`\mathrm{sin}^22\theta _{sun}`$. Various fits have been performed for two , three  and even four  neutrino oscillation scenarios. The most recent which take into account the full data (rates, energy spectrum, day-night asymmetry in the case of the MSW solution, and the seasonal variation in the case of the VO solution) give a best fit value of $`\mathrm{sin}^22\theta _{sun}<1`$. The $`95\%c.l.`$ bounds for the MSW solution are given in Table II. The most important conclusions are that in neither case do $`\mathrm{sin}^22\theta _{sun}`$ exceed the bound given by condition Eq. 8.

For the Vacuum Oscillation solution best fit value we take $`\mathrm{sin}^22\theta _{sun}=.93`$ . The above values have been used to draw Fig. 2 and Fig. 3 for the case $`s_+^e=.01`$.

Discussion. The most stringent of the oscillation solutions when question is of the neutrino nature, is the SMA. Already with the present bound one must adopt the Dirac description in the $`B_4`$ scheme. With negative results of GENIUS the same will hold true of the $`B_3`$ scheme. Currently a bound of $`m_\nu >.22eV`$ can be taken as a delimiter of the two neutrino natures for the $`A_3,A_4`$ and $`B_3`$ schemes.

At the $`95\%c.l.`$ the remaining solutions LMA and LOW define the neutrino nature if $`m_\nu >1.5eV`$, with the present bound and the $`A_3,A_4,B_3`$ schemes and $`m_\nu >1.1eV`$ in $`B_4`$. With the GENIUS I bound this can be shifted down to $`m_\nu <.16eV`$ ($`A_3`$), $`m_\nu <.14eV`$ ($`B_3`$) and $`m_\nu <.22eV`$ ($`A_4`$). Finally GENIUS II would give $`m_\nu <.05eV`$ for the $`A_3`$ and $`m_\nu <.12eV`$ for the $`A_4`$ scheme, but would completely exclude a Majorana neutrino in the $`B_3`$ scheme.

We do not have a $`95\%c.l.`$ bound for the VO solution. However, with the assumption that here also $`\mathrm{sin}^22\theta _{sun}<.98`$ the same conclusions would hold as above.

If one takes the best fit values for the different solutions, all of the bounds can be improved.

An interesting consideration is to take into account the current best fit for the Hot Dark Matter neutrino constituent. One then has $`_{light}m_\nu 2eV`$ . The SMA solution rules out in this case Majorana neutrinos in the $`A_3`$ and $`B_3`$ schemes. With a negative result from GENIUS all solutions share this property in the $`A_3,B_3`$ and $`B_4`$ schemes. Obviously, if we take the previous result $`_{light}m_\nu 45eV`$ , then the constraint is even stronger as already now most of the schemes are ruled out and GENIUS would close the situtation.

Final considerations. Several points must be stressed. First, the bound on $`|m_\nu |`$ depends strongly on the precise determination of nuclear matrix elements. It is acknowledged that they suffer from an imprecision at the level of a factor of 3 . Second, the confidence level that we assumed, $`95\%`$ corresponds to $`2\sigma `$. At $`3\sigma `$ which is $`99\%c.l.`$, the $`\mathrm{sin}^22\theta _{sun}`$ value of one is allowed. This unfortunately ruins the validity of condition Eq. 8 with the same effect on our predictions. Third, the figures have been drawn for $`\delta m_{LSND}^2=1eV^2`$. The experimentally allowed range is $`\delta m_{LSND}^2=.22eV^2`$. Adoption of some other difference of mass squared would require a rescaling of the results for the $`A_4`$ and $`B_4`$ neutrino schemes.

Conclusions. The nature of neutrinos that seemed so evasive is closer and closer to the reach of experiments. Combining oscillation data, tritium beta decay and neutrinoless double beta decay is a powerful method of constraining the parameter space in which Majorana neutrinos are allowed. In this paper we have derived bounds on the masses in several schemes and shown which schemes are or will be excluded out of consideration if neutrinos are to be truly neutral particles. We still have to wait for MiniBoone or SNO, to decide on the existence of the sterile neutrino, and for SNO, Borexino, Kamland to determine the solar neutrino solution and finally for more statistics in all the experiments to constrain more precisely the $`\mathrm{sin}^22\theta _{sun}`$.

Acknowledgments. This work was supported by the Polish Committee for Scientific Research under Grant No. 2P03B05418. J.G. would like to thank the Alexander von Humboldt-Stiftung for fellowship.


format: text

text:
# HARVESTING THERMAL FLUCTUATIONS: ACTIVATION PROCESS INDUCED BY A NONLINEAR CHAIN IN THERMAL EQUILIBRIUM

## 1 Introduction

The search for mechanisms that may induce the spontaneous localization of vibrational energy in molecular materials has surfaced in a variety of contexts where such localized energy may then trigger other events. These may include switching and other threshold phenomena, chemical reactions, local melting and other deformational effects, and even detonation. In the Kramers problem  a particle moving in a bistable potential is used as a model for a chemical process. The trajectory of the particle is associated with the reaction coordinate (RC). One well of the bistable potential represents the “reactant” state, the other the “product” state, and separating them is the “activation barrier.” The bistable potential is connected to a thermal environment, typically through fluctuating and dissipative terms, and every once in a while a large thermal fluctuation causes the particle to surmount the barrier and move from one well to the other. The average rate of occurrence of these events is associated with the reaction rate. This mesoscopic Langevin-type of approach admits of an underlying microscopic description of the thermal environment and its coupling to the bistable system. For instance, the usual Langevin equation with an instantaneous dissipation and Gaussian $`\delta `$-correlated fluctuations can be derived from a picture in which the system is harmonically coupled to an infinite number of harmonic oscillators with a uniform spectrum. A generalized Langevin picture involving dissipative memory terms and correlated fluctuations is associated with a more complex spectrum . It is clear, and has become a topic of considerable interest, that the nature of the environment and its coupling to the bistable system profoundly influence the transition rate.

A different but related set of problems that has attracted intense interest in recent years concerns the spontaneous localization of vibrational energy in periodic nonlinear arrays. The pioneering work of Fermi, Pasta and Ulam  demonstrated that a periodic lattice of coupled nonlinear oscillators is not ergodic, and that energy in such a lattice may never be distributed uniformly. A great deal of work has since followed in an attempt to understand how energy is distributed in discrete nonlinear systems . The existence of solitons and more generally of breathers and other energy-focusing mechanisms, and the stationarity or periodic recurrence or even slow relaxation of such spatially localized excitations, are viewed as nonlinear phenomena with important consequences in many physical systems . The search for localization mechanisms that are robust even when the arrays are in a thermal environment  has, on the one hand, narrowed the problem (because some localization mechanisms are fragile against thermal fluctuations) but on the other hand broadened it (because new entropy-driven localization mechanisms become possible). Thermal effects may be particularly important in biophysical and biochemical applications at the molecular level .

The interest in the distribution and motion of energy in periodic arrays arises in part because localized energy in these systems may be mobile, in contrast with systems where energy localization occurs through disorder. Localized energy that moves with little or no dispersion may appear at one location on an array and may then be able to move to another where it can be used in a subsequent process. Traditional harmonic models suffer from the fact that dispersion thwarts such a mechanism for energy transfer. There has been a surge of recent activity in an attempt to understand the thermal conductivity of nonlinear chains

The connection between the study of perfect nonlinear arrays and the Kramers problem arises because such arrays may themselves serve as models for a heat bath for other systems connected to them . Albeit in different contexts, “perfect” arrays serving as energy storage and transfer assemblies for chemical or photochemical processes are not uncommon , and literature on the subject goes back for two decades . We thus consider the following variant of the Kramers problem: a bistable system connected to a nonlinear chain, which is in turn connected to a heat bath in the usual Langevin manner (see Fig. 1). The bistable system is only connected to the environment through its embedding in the nonlinear chain, and therefore the ability of the chain to spontaneously localize thermal energy and perhaps to transport it to the location of the bistable system can profoundly affect the transition rate. We investigate the behavior of this model for different types of anharmonic chains and thereby establish the important role of the nature of the environment on these chemical model systems.

In Section 2 we discuss the energy landscape typical of various nonlinear chains in thermal equilibrium. In Section 3 the variant of the Kramers problem wherein a bistable system is connected to each of the different chains is presented. Section 4 details our results for the transition statistics in the bistable system. We compare and contrast the transition statistics in the different chains and compare them to those found in the standard (Markovian) Kramers and generalized Kramers problems. We conclude with a summary and some notes on future directions in Section 5.

## 2 Nonlinear Chains

The simplest nonlinear periodic arrays consist of masses connected by springs that may be harmonic or anharmonic. The masses may also experience a local harmonic or anharmonic potential. In a recent paper we presented a detailed view of the thermal landscape of arrays with local hard (the “$`\varphi ^4`$ model”), harmonic, or soft potentials and harmonic interactions . Here we present the complementary analysis (more interesting, it turns out, in the context of the Kramers problem) of the thermal landscape of masses connected by anharmonic springs (with no local potentials).

The spontaneous localization of energy in any system in thermal equilibrium is simply a reflection of the thermal fluctuations described by statistical mechanics and is unrelated to system dynamics. On the other hand, the way in which these fluctuations dissipate and/or move and disperse, that is, the temporal evolution of thermal fluctuations, is dictated by the system dynamics and, in particular, by the channels connecting the chain to the thermal environment (dissipation) and the masses to one another (intermolecular interactions).

We pose the following questions: 1) How is the energy distributed in an equilibrium nonlinear chain at any given instant of time, and how does this distribution depend on the anharmonicity? Can one talk about spontaneous energy localization in thermal equilibrium, and, if so, what are the mechanisms that lead to it? 2) How do local energy fluctuations in such an equilibrium array relax in a given oscillator? Are there circumstances in the equilibrium system wherein a given oscillator remains at a high level of excitation for a long time? 3) Can local high-energy fluctuations move in some nondispersive fashion along the array? Can an array in thermal equilibrium transmit long-lived high-energy fluctuations from one region of the array to another with little dispersion?

In our earlier work  we showed that in harmonically coupled nonlinear chains (“diagonal anharmonicity”) in thermal equilibrium, high-energy fluctuation mobility does not occur beyond that which is observed in a harmonic chain. Herein we show that the situation might be quite different if there is “nondiagonal anharmonicity,” that is, if the interoscillator interactions are anharmonic.

Our model consists of a one-dimensional array of $`N`$ unit-mass sites, each connected by a potential $`V(x_nx_{n\pm 1})`$ to its nearest neighbors that may be harmonic or anharmonic:

$$H=\underset{n=1}{\overset{N}{}}\frac{p_n^2}{2}+\underset{n=1}{\overset{N}{}}V(x_nx_{n1}).$$
(1)
We assume periodic boundary conditions and consider three prototype potentials:

$`V_h(x)`$ $`={\displaystyle \frac{1}{2}}kx^2+{\displaystyle \frac{1}{4}}k^{}x^4`$ $`\mathrm{hard},`$ (2)
$`V_0(x)`$ $`={\displaystyle \frac{1}{2}}kx^2`$ $`\mathrm{harmonic},`$ (3)
$`V_s(x)`$ $`={\displaystyle \frac{k}{k^{}}}\left[|x|{\displaystyle \frac{1}{k^{}}}\mathrm{ln}(1+k^{}|x|)\right]`$ $`\mathrm{soft}.`$ (4)
At small amplitudes the three potentials are harmonic with the same force constant $`k`$. The independent parameters $`k`$ and $`k^{}`$ allow control of the harmonic component and the degree of anharmonicity of the chain. Elsewhere  we have argued that the overarching characteristic of anharmonic oscillators is the dependence of frequency on energy. For a harmonic oscillator the frequency is $`\sqrt{k}`$ independent of energy; for a hard oscillator the frequency increases with energy, and for a soft oscillator the frequency decreases with energy.

The set of coupled stochastic equations of motion for the masses is that obtained from the Hamiltonian, Eq. (1), augmented by the usual Langevin prescription for coupling a system to a heat bath at temperature $`T`$:

$$\ddot{x}_n=\frac{}{x_n}\left[V(x_{n+1}x_n)+V(x_nx_{n1})\right]\gamma \dot{x}_n+\eta _n(t),$$
(5)
where a dot represents a derivative with respect to time. The $`\eta _n(t)`$ are mutually uncorrelated, zero-centered, Gaussian, $`\delta `$-correlated fluctuations that satisfy the fluctuation-dissipation relation $`\eta _n(t)\eta _j(t^{})=2\gamma k_BT\delta _{nj}\delta (tt^{})`$. The numerical integration of the stochastic equations for all our simulations is performed using the second order Heun’s method (which is equivalent to a second order Runge Kutta integration)  with time step $`\mathrm{\Delta }t=0.005`$. In each simulation the system is initially allowed to relax for enough iterations to insure thermal equilibrium, after which we take our “measurements.”

The equilibrium results to be presented here complement our observations, presented elsewhere, on the way in which these same chains propagate an energy pulse  as well as a sustained signal applied at a particular site .

A set of energy landscapes is shown in Fig. 2. Along the horizontal direction in each panel lies a thermalized chain of oscillators; the vertical upward progression shows the evolution of this equilibrium system with time. The gray scale represents the energy, with darker shading reflecting more energetic regions.

Several noteworthy features are evident in the figure. The energy fluctuations are greatest in the soft chain. This feature, seen earlier in chains with local anharmonic potentials , is a consequence of the effect that we have called entropic localization. In the soft chain not only are the thermal fluctuations greater at a given temperature, a result easily obtained from a simple virial analysis, but the free energy is minimized by a nonuniform distribution of energy that populates regions of phase space where the density of states is high. We have argued that this localization mechanism is robust against temperature increases – indeed, it becomes more effective with increasing temperature. A second distinctive feature of the soft chain is the persistence of the energy fluctuations: damping is not particularly effective for soft chains. The only other mechanism for removal of localized energy from a particular location is along the chain. This is clearly not an effective mechanism, a result that is in agreement with our analysis of the propagation of an externally applied pulse in the soft chain . The speed of propagation (in all chains) of a pulse of a given energy is essentially proportional to the average frequency associated with that energy, and in the soft chain this average frequency decreases with increasing energy . Although we do not see an obvious connection between these excitations and solitons at zero temperature (which are not entropic localization mechanisms) , there may be a closer connection with more generalized excitations such as breathers .

In the hard chain (Fermi-Pasta-Ulam chain) the total energy as well as the energy fluctuations are considerably smaller but quite mobile with little dispersion, as can easily be shown on the basis of a straightforward statistical mechanical analysis. In the hard chain the average frequency increases with energy and therefore more energetic pulses tend to travel more rapidly. We have also shown that the dispersion of energy in a hard chain is slow  – this is seen here in the integrity of the spontaneously localized pulses over a much longer time than in the harmonic chain.

An important question of course concerns the parameter regimes where the differences illustrated in Fig. 2 are observed. We have chosen potential parameters that insure clear distinctions in the displacement amplitudes associated with the three potentials at the chosen temperature. The only restriction is that the temperature not be “too low,” that is, we avoid the region where all three potentials are essentially harmonic. We have chosen very low damping for the illustration. The soft energy landscape is far less sensitive to damping than the hard array. An increase in damping would readily slow down the motion of the high-energy fluctuations and would shorten their lifetime. Further, while the speed of the energy fluctuation pulses is sensitive to the potential parameters, their lifetime and dispersion properties are less so (as long as one is in the highly anharmonic regime). On the other hand, the persistence of the fluctuations in the soft array is quite sensitive to the harmonic contribution to the potential. All else remaining fixed, the landscapes remain qualitatively similar as temperature increases: the fluctuations in the soft array become even stronger relative to the others, and the pulse speeds in the hard array become even higher.

Suppose now that a bistable “impurity” is embedded in each of these chains, as illustrated in Fig. 1. When sufficient energy reaches the impurity, a transition may occur from one well to the other. The statistical and dynamical properties of these transitions are not obvious, and are explored in the next section.

## 3 Kramers Problem and Statement of our Variant

### 3.1 Traditional Kramers Problem

The original Markovian Kramers problem  describes the reaction coordinate $`y`$ evolving in the bistable potential

$$V_b(y)=\frac{V_0}{4}(y^21)^2$$
(6)
according to the usual Langevin prescription

$$\ddot{y}=\frac{dV_b(y)}{dy}\gamma _b\dot{y}+\eta _b(t)$$
(7)
where $`\gamma _b`$ is the dissipation parameter (the subscript for “bistable” distinguishes this from the other dissipation parameters) and $`\eta _b(t)`$ represents Gaussian, zero-centered, $`\delta `$-correlated fluctuations that satisfy the fluctuation-dissipation relation appropriate for temperature $`T`$, $`\eta _b(t)\eta _b(t^{})=2\gamma _bk_BT\delta (tt^{})`$. The rate coefficient $`k_r`$ for transitions from one metastable well to the other is expressed as

$$k_r=\kappa k^{TST},$$
(8)
where $`k^{TST}`$ is the rate obtained from transition state theory for activated crossing, which in our units ($`V_0=1`$, frequencies $`\sqrt{2}`$ at the bottom of the two wells and unit frequency at the top of the barrier) is

$$k^{TST}=\frac{\sqrt{2}}{\pi }e^{1/4k_BT}.$$
(9)
This is the highest possible rate because it assumes no recrossings of the barrier when the particle moves from one well to the other. The “transmission coefficient” $`\kappa <1`$ captures the effects of recrossings. The dependence of $`\kappa `$ on the various parameters of the problem has been the subject of intense study over many years . Its dependence on $`\gamma _b`$ and temperature is exemplified in the simulations shown in Fig. 3. In particular, we note the occurrence of a maximum: as predicted by Kramers, the transmission coefficient at high friction (diffusion-limited regime) decays as $`\gamma _b^1`$ (and is independent of temperature); at low friction (energy-limited regime) Kramers predicted that $`\kappa `$ is proportional to $`\gamma _b/k_BT`$.

An important generalization of the original Kramers problem, the so-called Grote-Hynes problem , reformulates the model in terms of the generalized Langevin equation

$$\ddot{y}=\frac{dV_b(y)}{dy}_0^t𝑑t^{}\mathrm{\Gamma }(tt^{})\dot{y}(t^{})+\eta _b(t)$$
(10)
where $`\mathrm{\Gamma }(tt^{})`$ is a dissipative memory kernel and the fluctuation-dissipation relation is now generalized to $`\eta _b(t)\eta _b(t^{})=k_BT\mathrm{\Gamma }(tt^{})`$. The dissipative memory kernel reflects the dynamics of the thermal environment and is characterized by its own time scales. A frequent choice is an exponential, but other forms that have been used include a Gaussian and a decaying oscillatory memory kernel. The introduction of additional parameters of course changes the behavior of the transmission coefficient.

The main point here is to call attention to the fact that the transmission coefficient has the same value for two different values of the dissipation parameter $`\gamma _b`$, and that therefore one can not conclude whether the system is in one regime (diffusion-limited) or another (energy-limited) simply from the value of the transition rate. One needs to know the trend with changing dissipation parameter, and one requires further information about the dynamics underlying a given transition rate. Not surprisingly, these dynamics turn out to be entirely different in different regimes . The time dependence of the transmission coefficient is a direct reflection of the explicit trajectories of the particle as it transits from one well to the other. A number of investigators have looked at the time dependence of the transmission coefficient in the diffusion-limited  and energy-limited  regimes, and also at the effect of different types of memory kernels .

Of interest to us here are the different dynamical behaviors in the diffusion-limited regime and the energy-limited regime. In Fig. 4 we show two views of each of two typical trajectories of the reaction coordinate for the Markovian Kramers problem. The transmission coefficients associated with these two trajectories are not too different (see below), but

they correspond to different damping, placing them on opposite sides of the turnover in the $`\kappa `$ vs $`\gamma _b`$ curve. The trajectory in the first panel is in the diffusion-limited regime; the third panel shows an expanded view of a portion of this trajectory. The particle performs rather erratic motion within one well and once in a while it surmounts the barrier. When the particle surmounts the barrier it does not spend much time in the barrier region before being trapped again in one or the other well. The crossing trajectories thus tend to involve only one or a very small number of crossings/recrossings. The trajectory in the second panel, a portion of which is expanded in the fourth panel, is energy-limited. The particle performs a fairly periodic motion within one well. Barrier crossing events tend to retain the particle in the barrier region for several recrossings; a phase space analysis shows that the associated trajectories are rather smooth oscillations from one side to the other of the potential well above the barrier . Correlation functions associated with these trajectories are presented and discussed in Section 4.

### 3.2 Variant of the Kramers Problem

We would like to understand the way in which the very different thermal landscapes described in Section 2 affect the dynamics of a reaction coordinate evolving in a bistable “impurity” embedded in these environments. The connection of the bistable impurity to the thermal environment occurs only through its connection to the chain, that is, we set $`\gamma _b=0`$.

We need to specify how the bistable system interacts with the chain. We insert the impurity along the chain between sites $`i`$ and $`i+1`$ and connect it to each of these two sites (see Fig. 1). It is customary to choose a simple interaction potential with a harmonic dependence $`V_{int}(x,y)(xy)^2`$ for each chain site connected to the impurity. Here $`y`$ is the reaction coordinate and $`x`$ stands for the coordinate of the chain site connected to the impurity. However, this interaction tends to destabilize the bistability in that it causes the neighbors to pull the bistable particle toward the barrier rather than toward its natural metastable states. The interaction thus lowers the barrier of the bistable impurity. Since we do not want to “bias” the problem in this way, we have chosen an interaction that instead tends to favor the already metastable states:

$$V_{int}(x,y)=\frac{k_{int}}{2}\left(\frac{y^21}{2}x\right)^2.$$
(11)
Near the bistable minima (which are shifted by the interaction) the total potential for the reaction coordinate is still harmonic, and near the maximum at $`y=0`$ it is still parabolic. The barrier height is modulated by the motion of the neighbors (somewhat reminiscent of the barrier fluctuations in resonant activation problems). At large values of $`y`$ the interaction hardens the bistable potential.

The equations of evolution then have the following contributions. For a site in the chain not connected to the bistable impurity we have, as before, Eq. (5). For the bistable impurity

$$\ddot{y}=\frac{dV_b(y)}{dy}\frac{V_{int}(x_{i+1},y)}{y}\frac{V_{int}(x_i,y)}{y}.$$
(12)
For the site to the left of the bistable impurity

$$\ddot{x}_i=\frac{V(x_ix_{i1})}{x_i}\frac{V_{int}(x_i,y)}{x_i}\gamma \dot{x}_i+\eta _i(t)$$
(13)
and similarly for $`x_{i+1}`$. Comparisons are made for the same temperature, damping coefficients, and interaction parameter ($`k_{int}=0.1`$ throughout this work) varying only the nature of the chain.

Figure 5 shows trajectories of the bistable impurity embedded in each of the three chains. In the hard chain the trajectory is rather similar to that of a Markovian Kramers particle in the diffusion-limited regime, while in the soft chain it is closer to that of the energy-limited regime. This is a direct reflection of the behavior seen in Fig. 2, that is, of the fact that in the hard chain independent thermal fluctuations created elsewhere along the chain have a good chance of reaching the bistable impurity, causing erratic motion. An occasional large fluctuation causes a transition over the barrier, usually unaccompanied by recrossings: the same energy mobility that brings independent fluctuations to the impurity also makes it easy for the impurity to then lose a particular energy fluctuation back to the chain. In the soft chain, on the other hand, the particle performs fairly periodic motion within one well. Only fluctuations in the sites immediately adjacent to the impurity can excite the impurity; fluctuations originating elsewhere do not travel to the impurity. Strong fluctuations are therefore rarer but more energetic and more persistent, so transition events occur less often. However, once such a fluctuation occurs it tends to remain in that region for a long time; the reaction coordinate therefore recrosses the barrier a large number of times until it eventually loses the excess energy and is trapped again in one of the wells.

A second set of trajectories associated with the same bistable impurity in the same three chains at the same temperature but with a (10-fold) higher dissipation parameter is shown in Fig. 6. Not surprisingly, the trajectories are now more similar to one another, but nevertheless there are still important and revealing differences that will be made evident in our discussion in the next section. Furthermore, a comparison of the two sets will allow important observations concerning the trends associated with increased damping.

In the next section we provide a quantitative characterization of the differences in the trajectories and a comparison of these results with those of the traditional Kramers problem.

## 4 Results for Transition Rates

A useful description of the bistable system in different regimes is provided by the normalized correlation function

$$C(\tau )\frac{y(t+\tau )y(t)}{y^2(t)}$$
(14)
where the brackets indicate an average over $`t`$. Since $`y(t)=0`$, this correlation function decays to zero. When the thermal environment strongly and rapidly changes the particle momentum, the trajectory is erratic and the correlation function decays monotonically and exponentially. The decay time is a measure of the mean time between crossing events from one well to the other, and its inverse can be identified with the transition rate $`k_r`$. If on the other hand the effects of the thermal environment are weak, then the trajectory is determined mainly by the deterministic potential and remains correlated over much longer periods of time.

The correlation functions for the Markovian Kramers trajectories of Fig. 4 are shown in Fig. 7. In the high-dissipation regime the correlation function is monotonic and decays exponentially over essentially all times. This is a reflection of the essentially random motion within each well and between wells (the correlation functions for portions of the trajectory entirely within one well are also monotonically decreasing, albeit not to zero). The slope of the high-$`\gamma _b`$ curve in the right panel leads to a mean time between crossings of $`\tau _c250`$.

The oscillations in the low-$`\gamma _b`$ correlation function reflect mainly the systematic periodic motion of the particle within each well, i. e., of the portions of the trajectory that evolve for a long time near $`y=1`$ or near $`y=1`$. The period of these oscillations for the parameters used here is $`t_{bottom}=\sqrt{2}\pi `$, and this is very nearly the period of the oscillations in the figure. Crossing events from one well to the other are mostly separated by long times and are essentially independent (however, see further discussion below). Hence the logarithmic rendition in the right panel gives a straight line. Its slope leads to a mean time between crossing events of $`\tau _c450`$.

The Kramers correlation functions serve as a point of reference for an interpretation of the correlation functions associated with our variant of the Kramers problem. These are shown in a number of figures starting with Fig. 8, which shows the correlations functions associated with the trajectories in Figs. 5 and 6. We stress that in each panel the $`\gamma `$ and $`k_BT`$ are the same in all cases, as is the coupling of the chain to the bistable system; only the nature of the chain has changed. The first panel shows a correlation function for the harmonic chain that is oscillatory at early times, and quite similar to the energy-limited Markovian Kramers case (see also the corresponding trajectories in Figs. 5 and 4). We conjecture that the harmonic chain provides a thermal environment comparable to the low-damping Markovian Kramers environment. The correlation function associated with the hard chain is similar to the behavior at higher damping in the Kramers case. The correlation function associated with the soft chain also decays in an oscillatory fashion, but in a more complex way than in the energy-limited Markovian Kramers case. The alternation in the amplitudes is a consequence of the presence of sustained bursts of energy that cause a finite fraction of the trajectory to occur above the barrier, leading to many correlated recrossing events. The particle oscillates above the barrier for intervals much longer than in the Markovian Kramers trajectory. The typical oscillation period above the barrier is about twice as long as $`t_{bottom}`$ in our example (detailed discussions of these times can be found in our earlier work ). This effect is already slightly visible in the low-$`\gamma _b`$ Kramers correlation function in Fig. 7, but it is much stronger in the soft chain. To reproduce this behavior in the Kramers model it is necessary to consider the generalized Kramers model with a memory friction: there is clearly an additional memory effect in the soft chain that allows the energy to remain trapped in the region of the bistable impurity for a long time. This is in accord with the notion that transitions in the soft chain are caused by local thermal fluctuations in the nearest neighbors of the impurity. The impurity may periodically exchange energy with these neighbors before the fluctuation eventually dissipates away, and this causes repeated recrossings. This in turn leads to the conclusion that the memory friction in the Kramers model needed to reproduce the soft chain environment would most likely be oscillatory .

The periods of oscillations in the harmonic and soft chains are somewhat different from those of the Kramers curve in Fig. 7 and from each other. This is due to differences in the effective potentials.

The eventual decay of the correlation functions for all three chains is exponential. We find for the times between transition events (single or bursts as appropriate) $`\tau _c101`$ for the hard chain, $`\tau _c453`$ for the harmonic chain, and $`\tau _c1540`$ for the soft chain.

A similar set of correlation functions associated with the higher-damping-parameter trajectories of Fig. 6 is shown in the right panel of Fig. 8. The dynamics of the bistable system in the hard chain with increasing damping does not change in character, whereas the oscillations in the harmonic and soft chains become less pronounced as these systems move toward the diffusion-limited regime. First we note that all the curves become steeper, which translates to a shorter mean time between crossing events and therefore a higher transition rate for all three chains. The specific values we obtain are $`\tau _c72`$ for the hard chain, $`\tau _c170`$ for the harmonic chain, and $`\tau _c215`$ for the soft chain. The decrease in time between crossing events is most pronounced for the soft chain. This is consistent with the notion that the soft chain is in the energy-limited regime where small increases in effective damping cause the greatest increases in the transition rate (see Fig. 3). The harmonic chain lies closer in this sense to the turnover region, and the hard chain even closer yet. The second point is that this apparent trend for the hard chain indicates that it, too, lies on the low-damping side of the turnover in spite of the diffusion-limited aspects of its dynamics. This is the reason for the very small oscillations visible at the earliest times in the hard chain correlation functions. It is apparent that neither the trajectory itself nor even the shape of the correlation function at one value of the damping provides unequivocal information to determine which side of the turnover regime one is on; it is necessary to investigate the trend.

It is interesting to investigate whether our variants of the thermal environment can actually be “pushed” across the turnover point by increasing the damping on the chain. For this purpose we present a series of correlation functions for each of the chains for different values of the damping parameter. The other parameters, including the temperature, remain fixed and equal to the values given earlier.

Figure 9 shows results for the soft and harmonic chains. With increasing damping the early-time oscillations in the correlation function in the soft chain first lose some of the “alternation” feature typical of a long oscillatory dissipative memory kernel and eventually the correlation function loses its oscillatory character altogether. The crossing rate continues to increase as damping increases, so throughout this series one is still on the low damping side of the turnover. The right panel of Fig. 9 shows the correlation functions for the impurity in the harmonic chain. The trend is similar to that of the soft chain but, in all respects, indicative of the fact that the harmonic environment is closer to the turnover region than the soft environment. Thus the oscillations disappear sooner, and the increase in the transition rate with increasing damping is smaller.

Perhaps the most interesting features are seen in Fig. 10. Here we clearly see the very small short-time oscillations, which disappear as damping increases. The transition rate is quite insensitive to damping in the range $`\gamma =0.005`$$`0.1`$ shown in the figure (the line for $`\gamma =0.2`$, not explicitly shown, also falls in the same regime). The turnover value must therefore be within this range. To ascertain if this is so, we also exhibit the correlation function for a considerably larger value of the dissipation parameter, $`\gamma =1.0`$. The slower decay for $`\gamma =1`$ is clear in both panels.

## 5 Conclusions

There has been a dearth of information on the effects on the activation process of nonlinearities in the environment. We have taken an approach here that goes part way, much in the tradition of modeling efforts for a variety of systems interacting with a complex environment: the “immediate surroundings” of the reaction coordinate are described microscopically, while the interaction of this immediate environment with other degrees of freedom is handled phenomenologically.

We find that the dynamics of the activation process in some parameter regimes are profoundly affected by the nature of the chain. If the damping parameter connecting the chain to the heat bath is sufficiently low, a soft chain provides an environment very similar to that of the Grote Hynes model with an oscillatory memory kernel in the energy-limited regime , while a hard chain provides an environment akin to that of the Kramers model in the diffusion-limited regime . This in turn means that in such parameter regimes the hard chain is a more effective mediator of the activation process than is the soft chain.

A number of interesting questions concerning these systems are currently under investigation. One concerns the influence of boundary conditions on the behavior that we have described . A second problem concerns the effect on the reaction coordinate of a pulse or a sustained signal applied somewhere else along the chain. We have showed that the propagation of such a pulse or signal is strongly affected by the nature of the chain , and we expect these differences in turn to affect the response of a bistable impurity to these excitations. Such models are interesting in the context of physical or biophysical situations wherein energy is released at some location (provided perhaps by a chemical reaction or an absorption process at that location), which must then move to another location (that of the bistable impurity) to effect some further chemical process (represented by the activation process). The usual linear chain models are plagued by the excessive dispersion that would make such transmission inefficient. Nonlinearities in the environment may provide the necessary mobility with little attendant dispersion, thus greatly increasing the efficiency of such a process.

## Acknowledgments

R. R. gratefully acknowledges the support of this research by the Ministerio de Educación y Cultura through Postdoctoral Grant No. PF-98-46573147. A. S. acknowledges sabbatical support from DGAPA-UNAM. This work was supported in part by the Engineering Research Program of the Office of Basic Energy Sciences at the U.S. Department of Energy under Grant No. DE-FG03-86ER13606, and in part by the Comisión Interministerial de Ciencia y Tecnología (Spain) Project No. DGICYT PB96-0241.


text:
# The canonical transformations of the dynamical multiparameter systems as recurrence relations for the models on the grating

## I Introduction

Many equations of the theoretical physics especially of quantum mechanics are linear. In this connection there is a possibility of constructing solutions of these equations searching the canonical transformations of the linear dynamical systems (but this canonical transformations are different from those of quantum theory). These transformations can be formulated in terms of the creation and annihilation operators for some “quantum numbers” of the corresponding quantities that conserve their values. The indicated possibility was considered in , where the method of constructing the recurrence relation for searching for the eigenvalues and eigenfunctions of linear operators is proposed. This method is restricted by the case of one-dimensional systems described by the Lagrangian containing the spectral parameters in a potential part only.

In the present paper the generalization of this method for the case of the multidimensional systems and systems with the Lagrangian containing spectral parameters in a kinetic part too (for example, for the radial equation in hydrogen-like atom theory) is constructed.

The physical meaning of the considered canonical transformations for the classical theory differs from that for the quantum theory. The general point is that these transformations can be considered as translations on the multidimensional grating in the phase space. In the classical theory the parameters of the grating’s knots define the different configurations of considered class of dynamical systems (for example, the set of mechanical oscillators with different masses and spring rigidities). Conditions of the invariance under the canonical transformations along the grating define the admissible grating parameters and, therefore, solutions and a spectrum of oscillations.

In quantum theory the parameters of the knots give the possible states of the same quantum system. In its turn the translation along the grating is defined by the creation and annihilation operators stipulating the transition from one state to another.

The canonical recurrence transformations obtained in this way are the analogies of the known Beclund’s  transformations for finite-dimensional dynamical systems. Note that they do not coincide by their form with the standard recurrence relations cited in handbooks. The difference between them is explained by the way of their construction. The main here is that the canonical recurrence relations conserve the symplectic structure of the corresponding dynamical systems in contrast with the standard ones. That is why from the standpoint of the theory of dynamical systems the canonical recurrence relations have a certain advantage. This relations can be reduced to the standard ones by the corresponding normalization of eigenfunctions.

## II A canonical map

Let us consider the sets $`𝒟=\{q,p,H\}`$ and $`\stackrel{~}{𝒟}=\{\stackrel{~}{q},\stackrel{~}{p},\stackrel{~}{H}\}`$ of the $`s`$-dimensional dynamical systems with the generalized coordinates $`q=q(t)\{q^1(t),..,q^s(t)\}`$, $`\stackrel{~}{q}=\stackrel{~}{q}(t)\{\stackrel{~}{q}^1(t),\mathrm{},\stackrel{~}{q}^s(t)\}`$, momenta $`p=p(t)\{p_1(t),\mathrm{},p_s(t)\}`$, $`\stackrel{~}{p}=\stackrel{~}{p}(t)\{\stackrel{~}{p}_1(t),\mathrm{},\stackrel{~}{p}_s(t)\}`$ and Hamiltonians $`H=H(p,q,t)`$, $`\stackrel{~}{H}=\stackrel{~}{H}(\stackrel{~}{p},\stackrel{~}{q},t)`$.

We shall call the mapping $`:𝒟\stackrel{~}{𝒟}`$ the canonical one, if it is defined by the set $`=\{F\}`$ of generating functions $`F=F(q,\stackrel{~}{q},t)`$, so that for an arbitrary system $`(q,p,H)𝒟`$ and for some function $`F`$ there is the system $`(\stackrel{~}{q},\stackrel{~}{p},\stackrel{~}{H})\stackrel{~}{𝒟}`$, for which the relations

$$p_i=\frac{F}{q^i},\stackrel{~}{p}_i=\frac{F}{\stackrel{~}{q}^i},\stackrel{~}{H}=H+\frac{F}{t}(i,j,k,\mathrm{}=1,2,\mathrm{},s),$$
(1)
are realized.

If $`𝒟=\stackrel{~}{𝒟}`$, we shall say that given set of dynamical systems $`𝒟`$ is invariant under the canonical transformation $`F`$.

Let $`𝒟`$ be a parametrized class of dynamical systems described by the Lagrangian

$$L=L(q,\dot{q},t,\lambda ,\mu _a)=\frac{1}{2}\left[g_{ij}\dot{q}^i\dot{q}^j(U_{ij}\lambda g_{ij})q^iq^j\right],$$
(2)
where $`\dot{q}=dq/dt`$, $`g_{ij}=g_{ij}(\mu _a,t)`$, $`U_{ij}=U_{ij}(\mu _a,t)`$ are some functions, $`detg_{ij}0`$, $`\{\lambda ,\mu _a\}D`$, $`D`$ is some set of parameters, $`a=1,\mathrm{},d`$. The ordinary summation convention works here and on.

If we write the Lagrange-Euler equations for the system (2) in the form

$$\widehat{P}_k^lq^kg^{li}\left(\frac{d}{dt}g_{ik}\frac{d}{dt}+U_{ik}\right)q^k=\lambda q^l$$
(3)
we shall make a conclusion that $`\lambda `$ is a spectral parameter of the eigenvalue problem of the $`d`$-parametric set of the linear operators $`\widehat{P}_k^l(\mu _a,t)`$. For the Hamiltonians of the set (2) we have

$$H(p,q,t,\lambda ,\mu _a)=\frac{1}{2}g^{ij}p_ip_j\frac{1}{2}(U_{ij}\lambda g_{ij})q^iq^j,$$
(4)
where

$$p_i=g_{ij}\dot{q}^j,(g^{ik}g_{kj}=\delta _j^i).$$
(5)
The condition of the invariance of the set $`𝒟`$ of dynamical systems under the canonical transformations $``$ is

$$\stackrel{~}{H}(\stackrel{~}{p},\stackrel{~}{q},t,\stackrel{~}{\lambda },\stackrel{~}{\mu }_a)=H(\stackrel{~}{p},\stackrel{~}{q},t,\stackrel{~}{\lambda },\stackrel{~}{\mu }_a),\{\stackrel{~}{\lambda },\stackrel{~}{\mu }_a\}D.$$
(6)
Under sequential actions of the series of the canonical transformations a sequence of the systems of the type (2) and the parameters $`\{\lambda _m,\mu _{n_a}|m,n_a=1,2,\mathrm{}\}D`$ appears. Introduce the collective indices $`K=(m,n_a)`$, $`\stackrel{~}{K}=(\stackrel{~}{m},\stackrel{~}{n}_a)`$ and designate $`(\lambda _m,\mu _{n_a})\sigma _K`$, $`(\lambda _{\stackrel{~}{m}},\mu _{\stackrel{~}{n_a}})\sigma _{\stackrel{~}{K}}`$. Let the corresponding coordinates and momenta be $`q_K,q_{\stackrel{~}{K}}`$ and $`p_K,p_{\stackrel{~}{K}}`$. Further, $`F_{K\stackrel{~}{K}}=F_{K\stackrel{~}{K}}(q_K,q_{\stackrel{~}{K}},t)`$ be the generating function of the canonical transformation $`F`$: $`𝒟𝒟`$, so that $`\{q_K,p_K,H(p_K,q_K,t,\sigma _K)\}`$ $``$ $`\{q_{\stackrel{~}{K}},p_{\stackrel{~}{K}},H(p_{\stackrel{~}{K}},q_{\stackrel{~}{K}},t,\sigma _{\stackrel{~}{K}})\}`$. Then the conditions (1), (6) give the following equation for the generating function

$$\frac{F_{K\stackrel{~}{K}}}{t}+H(\frac{F_{K\stackrel{~}{K}}}{q_K},q_K,\sigma _K,t)=H(\frac{F_{K\stackrel{~}{K}}}{q_{\stackrel{~}{K}}},q_{\stackrel{~}{K}},\sigma _{\stackrel{~}{K}},t).$$
(7)
Further we shall write the indices $`K`$ and $`\stackrel{~}{K}`$ for the generating functions only, and for an arbitrary physical variable $`f`$ we shall write $`f`$ or $`\stackrel{~}{f}`$ instead of $`f_K`$ or $`f_{\stackrel{~}{K}}`$ respectively. Then it is possible to rewrite the condition of the invariance (7) for the system (2) with the Hamiltonian (4) in the following form:

$$2\frac{F_{K\stackrel{~}{K}}}{t}+g^{ij}\frac{F_{K\stackrel{~}{K}}}{q^i}\frac{F_{K\stackrel{~}{K}}}{q^j}\stackrel{~}{g}^{ij}\frac{F_{K\stackrel{~}{K}}}{\stackrel{~}{q}^i}\frac{F_{K\stackrel{~}{K}}}{\stackrel{~}{q}^j}=(U_{ij}\lambda g_{ij})q^iq^j(\stackrel{~}{U}_{ij}\stackrel{~}{\lambda }\stackrel{~}{g}_{ij})\stackrel{~}{q}^i\stackrel{~}{q}^j.$$
(8)
## III  A recurrence relations as a canonical map

Let us consider the problem of searching for the solutions of Eq.(8). The canonical map, that converts the linear dynamical systems into linear ones, must be linear. That is why the generating function $`F_{K\stackrel{~}{K}}`$ must be a quadratic function of the generalized coordinates $`q,\stackrel{~}{q}`$. Therefore we search for $`F_{K\stackrel{~}{K}}`$ in the form

$$F_{K\stackrel{~}{K}}=\frac{1}{2}\left(2\gamma _{ij}q^i\stackrel{~}{q}^jb_{ij}q^iq^jc_{ij}\stackrel{~}{q}^i\stackrel{~}{q}^j\right),$$
(9)
where $`\gamma _{ij},b_{ij},c_{ij}`$ are some functions. Substituting Eq.(9) for Eq.(8) and equating the coefficients at $`q^iq^j`$, $`q^i\stackrel{~}{q}^j`$, $`\stackrel{~}{q}^i\stackrel{~}{q}^j`$, we obtain

$`2\dot{\gamma }_{rs}g^{ij}(\gamma _{ir}b_{js}+\gamma _{jr}b_{is})=\stackrel{~}{g}^{ij}(\gamma _{is}c_{jr}+\gamma _{js}b_{ir}),`$ (10)
$`\dot{b}_{rs}+g^{ij}b_{ir}b_{js}+\lambda g_{rs}U_{rs}=\stackrel{~}{g}^{ij}\gamma _{ir}\gamma _{js},`$ (11)
$$\dot{c}_{rs}+g^{ij}\gamma _{ir}\gamma _{js}=\stackrel{~}{g}^{ij}c_{ir}c_{js}\stackrel{~}{U}_{rs}+\stackrel{~}{\lambda }\stackrel{~}{g}_{rs}.$$
(12)
If we find one of the particular solutions $`\gamma _{ij},b_{ij},c_{ij}`$, then, using Eqs.(1), (5), (9), it will be easy to construct the binomial recurrence relations for the eigenfunctions $`q(t)`$ and $`\stackrel{~}{q}(t)`$ of the set of operators $`\widehat{P}_k^l`$

$$\left(b_{ij}+g_{ij}\frac{d}{dt}\right)q^j=\gamma _{ij}\stackrel{~}{q}^j,$$
(13)
$$\left(c_{ij}\stackrel{~}{g}_{ij}\frac{d}{dt}\right)\stackrel{~}{q}^j=\gamma _{ij}q^j.$$
(14)
The obtained relations, increasing and decreasing $`K`$, are the analogies of the well-known Beclund’s transformations  for finite-dimensional systems. They are also the generalization of the known relations of the factorization method , . It means, that Eqs.(13), (14) install the correlation between the canonical transformation method and factorization one.

The algebraic trinomial recurrence relations can be obtained considering the sum of the generating functions of the two sequential transformations

$$F_{K\stackrel{~}{K}\stackrel{~}{\stackrel{~}{K}}}=F_{K\stackrel{~}{K}}+F_{\stackrel{~}{K}\stackrel{~}{\stackrel{~}{K}}}.$$
(15)
According to Eq.(1) we have $`F_{K\stackrel{~}{K}\stackrel{~}{\stackrel{~}{K}}}/\stackrel{~}{q}=0`$, from which the linear algebraic recurrence relation

$$\gamma _{ij}q^j(c_{ij}+\stackrel{~}{b}_{ij})\stackrel{~}{q}^j+\stackrel{~}{\gamma }_{ij}\stackrel{~}{\stackrel{~}{q}}^j=0,$$
(16)
follows. Here $`\stackrel{~}{b}_{ij},\stackrel{~}{c}_{ij},\stackrel{~}{\gamma }_{ij}`$ are the coefficients of a quadratic form of the type (9) for the generating function $`F_{\stackrel{~}{K}\stackrel{~}{\stackrel{~}{K}}}`$.

The generating function $`F_{K\stackrel{~}{\stackrel{~}{K}}}`$ of the canonical transformation $`(q,p)(\stackrel{~}{\stackrel{~}{q}},\stackrel{~}{\stackrel{~}{p}})`$, which is a composition of the canonical transformations $`(q,p)(\stackrel{~}{q},\stackrel{~}{p})`$ and $`(\stackrel{~}{q},\stackrel{~}{p})(\stackrel{~}{\stackrel{~}{q}},\stackrel{~}{\stackrel{~}{p}})`$, can be constructed excluding the intermediate state $`\stackrel{~}{q}`$ from Eq.(15) by means of Eq.(16). It makes possible to simplify the procedure of searching for the generating function of the transformation $`K\stackrel{~}{K}`$ decomposing it into more simple steps of the calculating the elementary generating functions of the transformations: $`K=\{m,n_1,\mathrm{}\}\{m+1,n_1,\mathrm{}\}\{m+1,n_1+1,\mathrm{}\}\mathrm{}\stackrel{~}{K}`$, and after this to construct their composition.

From the standpoint of classical mechanics we consider here the canonical transformation theory of an infinite-dimensional system with the action

$$S=\frac{1}{2}𝑑t\underset{K}{}\left\{g_{ij}\dot{q}_K^i\dot{q}_K^j+(U_{ij}\lambda _mg_{ij})q_K^iq_K^j\right\},$$

describing a multicomponent model on a $`(d+1)`$-dimensional grating. The values of the parameters $`\sigma _K=\{\lambda _m,\mu _{n_1},\mathrm{},\mu _{n_d}\}`$ that define grating knots are determined from the condition of an invariance of the dynamical system under grating translations induced by the canonical transformations of the system.

One more interesting possibility connected with the method of the inverse scattering problem is shown below by the examples of more simple dynamical systems.

For the particular case of the set of $`d`$-component three-parametrical dynamical systems $`(s=1)`$ let us denote: $`g_{11}m(\nu ,t)`$, $`U_{11}m(\nu ,t)U(\mu ,t)`$, where $`\mu ,\nu `$ are parameters. Write the Lagrangian (2) in the form

$$L=\frac{m(\nu ,t)}{2}\left[\dot{q}^2+(U(\mu ,t)\lambda )q^2\right].$$
(17)
It is convenient to find the generating function $`F_{K\stackrel{~}{K}}`$ in such a form

$$F_{K\stackrel{~}{K}}=mS,S=\frac{1}{2a}(2\gamma q\stackrel{~}{q}bq^2c\stackrel{~}{q}^2),$$
(18)
where $`\gamma `$ is some unknown constant, $`a,b,c`$ are desired functions of $`t`$. The conditions (1) may be rewritten in the form

$$\dot{q}=\frac{S}{q},\dot{\stackrel{~}{q}}=\frac{m}{\stackrel{~}{m}}\frac{S}{\stackrel{~}{q}},\stackrel{~}{H}=H+\frac{}{t}(mS).$$
(19)
The recurrence relations

$$\left(b+a\frac{d}{dt}\right)q=\gamma \stackrel{~}{q},\left(ca\frac{m}{\stackrel{~}{m}}\frac{d}{dt}\right)\stackrel{~}{q}=\gamma q$$
(20)
and equation for the function $`S`$

$$\frac{2}{m}\frac{}{t}(mS)+\left(\frac{S}{q}\right)^2\frac{m}{\stackrel{~}{m}}\left(\frac{S}{t}\right)^2(\lambda U)q^2\frac{\stackrel{~}{m}}{m}(\stackrel{~}{\lambda }\stackrel{~}{U})\stackrel{~}{q}^2=0.$$
(21)
follow from Eq.(6), (19). Eqs. (10)–(12) for determining $`a,b,c,\gamma `$ acquire the form

$$\dot{a}=\frac{\dot{m}}{m}a+\frac{m}{\stackrel{~}{m}}cb,$$
(22)
$$\frac{m}{\stackrel{~}{m}}\dot{c}+\dot{b}=(\lambda \stackrel{~}{\lambda }+\stackrel{~}{U}U)a,$$
(23)
$$a\dot{b}=(\lambda U)a^2+\frac{m}{\stackrel{~}{m}}(bc\gamma ^2).$$
(24)
The system of equations (22)–(23) is linear with respect to the unknown functions $`a,b,c`$ and does not contain the constant $`\gamma `$, while the equation (24) is quadratic and contains $`\gamma `$. The procedure of the searching for the solutions of the system (22)–(24) consists in the finding a more simple particular solution of the homogeneous linear undetermined system (22)- (23) for which the quadratic equation (24) is satisfied at the some unknown parameters $`\{\stackrel{~}{\lambda },\stackrel{~}{\mu },\stackrel{~}{\nu }\}`$. Then the constant $`\gamma `$ will be determined.

The algebraic trinomial recurrence relations (without derivatives) and their compositions can be found by the way analogous to obtaining Eq.(16) (see  also). As a result we shall have

$$\frac{\gamma _{KK+1}}{a_K}q_K\left(\frac{c_{K+1K+1}}{a_K}+\frac{b_{K+1K+1}}{a_{K+1}}\right)q_{K+1}+\frac{\gamma _{K+1K+2}}{a_{K+1}}q_{K+2}=0,$$
(25)
where $`K`$ is the collective index of the totality of the parameters $`\{\lambda _m,\mu _n,\nu _l\}`$, i.e. if $`K=\{m,n,l\}`$ then $`K+1`$ means that one of the indices $`m,n,l`$ is increased by one. The coefficients in Eq.(25) correspond to that of the generating functions $`S_{KK+1}`$ and $`S_{K+1K+2}`$ of the type (18) of the transformations $`q_Kq_{K+1}`$ and $`q_{K+1}q_{K+2}`$.

Writing the similar equations for the sequence of the transformations $`q_{K+1}q_{K+2}q_{K+3}`$ and excluding the coordinate $`q_{K+2}`$ by means of Eq.(25), we obtain the two relations, which can be written in the matrix form

$$Q_{K+2}=B_KQ_K,$$
(26)
where

$$Q_K=\left(\begin{array}{c}q_K\\ q_{K+1}\end{array}\right),Q_{K+2}=\left(\begin{array}{c}q_{K+2}\\ q_{K+3}\end{array}\right),$$
(27)
and $`B_K`$ is a $`2\times 2`$ matrix, which is expressed through the coefficients of the quadratic forms $`\{S_{KK+1},\mathrm{},S_{K+2K+3}\}`$ of the above sequence of the transformations; its explicit form is not indicated here because of its inconvenience.

From the other hand the differential equation (20) can be rewritten in these terms as

$$\frac{dQ_K}{dt}=A_KQ_K,$$
(28)
where

$$A_K=\frac{1}{a_Km_{K+1}}\left(\begin{array}{cc}m_{K+1}b_{KK}& m_{K+1}\gamma _{KK+1}\\ m_K\gamma _{KK+1}& m_Kc_{K+1K+1}\end{array}\right).$$
(29)
The conditions of compatibility of Eqs.(26)–(27)

$$\frac{dB_K}{dt}+B_KA_KA_{K+2}B_K=0$$
(30)
are equations for the coefficients $`a_K,b_K,c_K,\gamma _K`$ of the sequence of the generating functions $`S_{KK+1}`$ and represent zero curvature conditions for the models on the grating . They are the conditions of an applicability of the inverse scattering problem method for integrable models. In our case they are realized as a consequence of the equations of the type (22)–(24).

Thus the regular way of the constructing zero curvature conditions, which can be generalized for nonlinear dynamical integrable systems, follows from our analysis.

## IV The confluent hypergeometric equation

Many equations of quantum mechanics under the separation of variables and after isolation of an angular part lead to the second order ordinary differential equations, which can be reduced to the hypergeometric (or confluent hypergeometric) equation. In this connection it is important to consider the canonical transformation theory of the equations of such a type in general case.

Let us consider construction of the canonical recurrence relations for the confluent hypergeometric equation

$$t\ddot{q}+(\beta t)\dot{q}\alpha q=0,$$
(31)
where $`\alpha ,\beta `$ are some parameters. It is the Lagrange-Euler equation for a dynamical system with the Lagrangian

$$L=\frac{1}{2}t^\beta e^t(\dot{q}^2+\frac{\alpha }{t}q^2).$$
(32)
Comparing Eq.(32) with Eq.(17), we can conclude that

$$\lambda =0,\nu =\beta ,\mu =\alpha ,U(\alpha ,t)=\frac{\alpha }{t},m(\beta ,t)=t^\beta e^t.$$

Therefore the system (22)–(24) has the form

$$\dot{a}=\frac{\dot{m}}{m}a+\frac{m}{\stackrel{~}{m}}cb,\frac{m}{\stackrel{~}{m}}\dot{c}+\dot{b}=\frac{\stackrel{~}{\alpha }\alpha }{t}a$$
(33)
$$a\dot{b}=\frac{\alpha }{t}a^2+\frac{m}{\stackrel{~}{m}}(bc\gamma ^2)$$
(34)
According to the above-mentioned, we shall find the generating functions of the elementary canonical transformations a) $`\{\alpha ,\beta \}\{\stackrel{~}{\alpha },\beta \}`$, b) $`\{\alpha ,\beta \}\{\alpha ,\stackrel{~}{\beta }\}`$ only. Besides note that Eqs.(33)–(34) can be written as equations with rational on $`t`$ coefficients. Therefore it is naturally to find the solution in a class of rational functions.

### A A canonical transformation $`\{\alpha ,\beta \}\{\stackrel{~}{\alpha },\beta \}`$

In this case $`\beta =\stackrel{~}{\beta }`$, $`m=\stackrel{~}{m}`$ and the system (33)–(34) acquires the form

$$\dot{a}=(\frac{\beta }{t}1)a+cb,\dot{c}+\dot{b}=\frac{\stackrel{~}{\alpha }\alpha }{t}a,$$
(35)
$$a\dot{b}=bc\gamma ^2a^2\frac{\alpha }{t}.$$
(36)
We shall obtain the simplest rational solution of the undetermined system (35) by assumption $`a=t`$. From (35) we shall find

$$b=\frac{1}{2}(\stackrel{~}{\alpha }\alpha 1)t+b_0,c=\frac{1}{2}(\stackrel{~}{\alpha }\alpha +1)t+b_0\beta +1,$$
(37)
where $`b_0`$ is the integration constant. Substituting Eq.(37) for Eq.(36) and equating the coefficients at the same powers of $`t`$, one obtains

$$b_0=\alpha ,\stackrel{~}{\alpha }=\alpha +1,\gamma =\sqrt{\alpha (\alpha \beta +1)}$$
(38)
We have restricted ourselves here by the positive value of the radical. As a result we have

$$F_{\alpha ,\beta \alpha +1,\beta }=\frac{1}{2}t^{\beta 1}e^t\left(2\sqrt{\alpha (\alpha \beta +1)}q\stackrel{~}{q}\alpha q^2(t+\alpha \beta +1)\stackrel{~}{q}^2\right)$$
(39)
for the generating function (18). The canonical recurrence transformations (20) have the form

$$\left(\alpha +t\frac{d}{dt}\right)q=\sqrt{\alpha (\alpha \beta +1)}\stackrel{~}{q},\left(t+\alpha \beta +1t\frac{d}{dt}\right)\stackrel{~}{q}=\sqrt{\alpha (\alpha \beta +1)}q.$$
(40)
To determine the connection with the standard recurrence relations and thus with the confluent hypergeometric function $`M(\alpha ,\beta ,t)`$, we shall make the following substitutions

$$q=y(\alpha )\sqrt{N(\alpha )},\stackrel{~}{q}=y(\alpha )\sqrt{N(\alpha +1)},y(\alpha )M(\alpha ,\beta ,t)$$
(41)
in Eq.(40). From the condition of coincidence of the obtained recurrence relations with the standard ones  (rather with their consequence)

$$\left(\alpha +t\frac{d}{dt}\right)y(\alpha )=\alpha y(\alpha +1),\left(\alpha \beta +t1t\frac{d}{dt}\right)y(\alpha +1)=(\alpha \beta +1)y(\alpha )$$
(42)
the functional equation for the normalization multiplier $`N(\alpha )`$

$$\alpha N(\alpha )=(\alpha \beta +1)N(\alpha +1)$$
(43)
follows. It can be solved by using the gamma-function’s property $`\mathrm{\Gamma }(z+1)=z\mathrm{\Gamma }(z)`$ from where the substitutions

$$\alpha =\frac{\mathrm{\Gamma }(\alpha +1)}{\mathrm{\Gamma }(\alpha )},\alpha \beta +1=\frac{\mathrm{\Gamma }(\alpha \beta +2)}{\mathrm{\Gamma }(\alpha \beta +1)}.$$
(44)
follow. Using these substitutions, equation (43) can be reduced to the form

$$\frac{\mathrm{\Gamma }(\alpha \beta +1)}{\mathrm{\Gamma }(\alpha )}N(\alpha )=\frac{\mathrm{\Gamma }(\alpha \beta +2)}{\mathrm{\Gamma }(\alpha +1)}N(\alpha +1)=\frac{\mathrm{\Gamma }(\alpha \beta +3)}{\mathrm{\Gamma }(\alpha +2)}N(\alpha +2)=\mathrm{}.$$
(45)
The additional equations arise owing to considering the next elementary steps by the parameter $`\alpha `$. From here it is easy to see that $`\mathrm{\Gamma }(\alpha \beta +1)N(\alpha )/\mathrm{\Gamma }(\alpha )=C_1=\mathrm{const}`$. Therefore

$$N(\alpha )=C_1\frac{\mathrm{\Gamma }(\alpha )}{\mathrm{\Gamma }(\alpha \beta +1)}.$$

Thus, we have obtained the “canonical” solution of the confluent hypergeometrical equation in the form

$$q(t)=C_1M(\alpha ,\beta ,t)\sqrt{\frac{\mathrm{\Gamma }(\alpha )}{\mathrm{\Gamma }(\alpha \beta +1)}}.$$
(46)
### B A canonical transformation $`\{\alpha ,\beta \}\{\alpha ,\stackrel{~}{\beta }\}`$

In this case $`\alpha =\stackrel{~}{\alpha }`$ and the system (33)–(34) can be rewritten in the form

$$t^{\beta \stackrel{~}{\beta }}cb=\dot{a}+(1\frac{\beta }{t})a,t^{\beta \stackrel{~}{\beta }}\dot{c}+\dot{b}=0,$$
(47)
$$a\dot{b}+\frac{\alpha }{t}a^2+t^{\beta \stackrel{~}{\beta }}(\gamma ^2bc)=0.$$
(48)
We shall obtain the simplest rational solution of the inmohomogeneous undetermined system (47) by taking $`a=1`$. Then a particular solution of the obtained nonhomogeneous system has the form

$$c=\frac{\beta }{\stackrel{~}{\beta }\beta +2}t^{\stackrel{~}{\beta }\beta 1},b=\frac{\stackrel{~}{\beta }\beta 1}{\stackrel{~}{\beta }\beta +2}\frac{\beta }{t}.$$
(49)
From the quadratic equation (48) we find

$$\stackrel{~}{\beta }=\beta +1,\gamma =\sqrt{\beta \alpha }.$$
(50)
For the generating function and canonical recurrence transformations we have

$$F_{\alpha ,\beta \alpha ,\beta +1}=\frac{1}{2}t^\beta e^t\left(2\sqrt{\beta \alpha }q\stackrel{~}{q}+q^2+\beta \stackrel{~}{q}^2\right),$$
(51)
$$\left(1+\frac{d}{dt}\right)q=\sqrt{\beta \alpha }\stackrel{~}{q},\left(\beta +t\frac{d}{dt}\right)\stackrel{~}{q}=\sqrt{\beta \alpha }q.$$
(52)
The recurrence relations increasing and decreasing $`\beta `$ that follow from the recurrence relations of the handbook

$$\beta \left(1+\frac{d}{dt}\right)y(\beta )=(\beta \alpha )y(\beta +1),\left(\beta +t\frac{d}{dt}\right)y(\beta +1)=\beta y(\beta ),$$
(53)
where $`y(\beta )M(\alpha ,\beta ,t)`$, can be obtained from Eq.(52) by the above procedure of the constructing of Eq. (42)–(46). As a result one has

$$q(t)=C_2\frac{\mathrm{\Gamma }(\beta )}{\sqrt{\mathrm{\Gamma }(\beta \alpha )}}M(\alpha ,\beta ,t).$$
(54)
## V The hypergeometric equation

As it was mentioned, some “radial equations” of the quantum theory are reduced to the hypergeometric-like equations. Therefore consider a dynamical system described by the hypergeometrical equation

$$t(1t)\ddot{q}((\alpha +\beta +1)t\zeta )\dot{q}\alpha \beta q=0,$$
(55)
where $`\alpha ,\beta ,\zeta `$ are some parameters, as another application of the developed theory. It is the Lagrange-Euler equation of a dynamical system with the Lagrangian

$$L=\frac{1}{2}t^\zeta (1t)^{\alpha +\beta \zeta +1}\left(\dot{q}^2+\frac{\alpha \beta }{t(1t)}q^2\right).$$
(56)
Comparing Eq.(56) with Eqs.(2), (17) we have: $`\lambda =0,\mu _1=\alpha ,\mu _2=\beta ,\mu _3=\zeta ,`$

$$U=\frac{\alpha \beta }{t(1t)},m=t^\zeta (1t)^{\alpha +\beta \zeta +1}.$$
(57)
The system (56) contains the three parameters $`\alpha ,\beta ,\zeta `$. Therefore it is necessary to find the generating functions of three canonical transformations. Due to the symmetry between the parameters $`\alpha `$ and $`\beta `$, it is sufficient to search for the generating functions of the transformations a) $`\{\alpha ,\beta ,\zeta \}\{\stackrel{~}{\alpha },\beta ,\zeta \}`$ and b) $`\{\alpha ,\beta ,\zeta \}\{\alpha ,\beta ,\stackrel{~}{\zeta }\}`$ only.

### A A canonical transformation $`\{\alpha ,\beta ,\zeta \}\{\stackrel{~}{\alpha },\beta ,\zeta \}`$

For the elementary transformation $`\{\alpha ,\beta ,\zeta \}\{\stackrel{~}{\alpha },\beta ,\zeta \}`$ the system (22)–(24) can be rewritten the form

$$(1t)^{\alpha \stackrel{~}{\alpha }}cb=\dot{a}\left(\frac{\zeta }{t}\frac{\alpha +\beta \zeta +1}{1t}\right)a,$$
(58)
$$(1t)^{\alpha \stackrel{~}{\alpha }}\dot{c}+\dot{b}=\frac{\stackrel{~}{\alpha }\alpha }{t(1t)}\beta a,$$
(59)
$$a\dot{b}=\frac{\alpha \beta }{t(1t)}a^2+(1t)^{\alpha \stackrel{~}{\alpha }}(bc\gamma ^2).$$
(60)
Similarly to the first case of the previous example we obtain the simplest rational solution of the homogeneous undetermined system (58)–(59) at $`a=t`$. Then the obtained inhomogeneous system has the following particular solution

$$b=\alpha \frac{\stackrel{~}{\alpha }\alpha 1}{\stackrel{~}{\alpha }\alpha 2}\frac{\alpha +\beta \zeta +1}{1t},c=\left(\beta +\frac{\alpha +\beta \zeta +1}{\stackrel{~}{\alpha }\alpha 2}\frac{1}{1t}\right)(1t)^{\stackrel{~}{\alpha }\alpha }.$$
(61)
In this case the quadratic equation (60) is satisfied at

$$\alpha =\stackrel{~}{\alpha }+1,\gamma =\sqrt{\alpha (\alpha \zeta +1)}.$$
(62)
So for the generating function we have

$$F_{\alpha ,\beta ,\zeta \alpha +1,\beta ,\zeta }=\frac{1}{2}t^{\zeta 1}(1t)^{\alpha +\beta \zeta +1}\left[2\sqrt{\alpha (\alpha \zeta +1)}q\stackrel{~}{q}\alpha q^2(\alpha \zeta +1+\beta t)\stackrel{~}{q}^2\right].$$
(63)
The corresponding recurrence transformations have the form

$`\left(\alpha +t{\displaystyle \frac{d}{dt}}\right)q`$ $`=`$ $`\sqrt{\alpha (\alpha \zeta +1)}\stackrel{~}{q},`$ (64)
$`\left(\alpha \zeta +1\beta tt(1t){\displaystyle \frac{d}{dt}}\right)\stackrel{~}{q}`$ $`=`$ $`\sqrt{\alpha (\alpha \zeta +1)}q.`$ (65)
The standard recurrence transformation

$`\left(\alpha +t{\displaystyle \frac{d}{dt}}\right)y(\alpha )`$ $`=`$ $`\alpha y(\alpha +1),`$ (66)
$`\left(\alpha \zeta +1+\beta tt(1t){\displaystyle \frac{d}{dt}}\right)\stackrel{~}{q}`$ $`=`$ $`\sqrt{\alpha (\alpha \zeta +1)}q,`$ (67)
where $`y(\alpha )M(\alpha ,\beta ,\zeta ,t)`$ — hypergeometric function, can be obtained by making such the substitution

$$q(t)=A_1\sqrt{\frac{\mathrm{\Gamma }(\alpha )}{\mathrm{\Gamma }(\alpha \zeta +1)}}M(\alpha ,\beta ,\zeta ,t)$$
(68)
in Eq.(64).

### B A canonical transformation $`\{\alpha ,\beta ,\zeta \}\{\alpha ,\beta ,\stackrel{~}{\zeta }\}`$

In this case the system (22)–(24) acquires the form

$`t^{\zeta \stackrel{~}{\zeta }}(1t)^{\stackrel{~}{\zeta }\zeta }cb`$ $`=`$ $`\dot{a}\left({\displaystyle \frac{\zeta }{t}}{\displaystyle \frac{\alpha +\beta \zeta +1}{1t}}\right)a,`$ (69)
$`t^{\zeta \stackrel{~}{\zeta }}(1t)^{\stackrel{~}{\zeta }\zeta }\dot{c}+\dot{b}`$ $`=`$ $`0,`$ (70)
$$a\dot{b}=\frac{\alpha \beta }{t(1t)}a^2+t^{\zeta \stackrel{~}{\zeta }}(1t)^{\stackrel{~}{\zeta }\zeta }(b\gamma ^2).$$
(71)
We obtain the simplest solution of the system (69) at $`a=1t`$. Then a particular solution (69) acquires the form

$$b=\frac{\zeta }{\stackrel{~}{\zeta }\zeta 2}\frac{1t}{t}+\frac{\zeta }{t}\alpha \beta ,c=\frac{\zeta }{\stackrel{~}{\zeta }\zeta 2}t^{\stackrel{~}{\zeta }\zeta 1}(1t)^{\zeta \stackrel{~}{\zeta }+1}.$$
(72)
The quadratic equation (71) leads to the following constraints

$$\stackrel{~}{\zeta }=\zeta +1,\gamma =\sqrt{(\zeta \alpha )(\beta \zeta )}.$$
(73)
Therefore for the generating function and the recurrence relations we find

$$F_{\alpha ,\beta ,\zeta \alpha ,\beta ,\zeta +1}=\frac{1}{2}t^\zeta (1t)^{\alpha +\beta \zeta }\left(2\sqrt{(\zeta \alpha )(\beta \zeta )}q\stackrel{~}{q}+(\alpha +\beta \zeta )q^2+\zeta \stackrel{~}{q}^2\right),$$
(74)
$`\left(\zeta \alpha \beta +(1t){\displaystyle \frac{d}{dt}}\right)q`$ $`=`$ $`\sqrt{(\zeta \alpha )(\beta \zeta )}\stackrel{~}{q},`$ (75)
$`\left(\zeta +t{\displaystyle \frac{d}{dt}}\right)\stackrel{~}{q}`$ $`=`$ $`\sqrt{(\zeta \alpha )(\beta \zeta )}q.`$ (76)
Making the substitution

$$q(t)=A_2\sqrt{\frac{\mathrm{\Gamma }(\zeta \alpha )\mathrm{\Gamma }(\zeta \beta )}{\mathrm{\Gamma }(\zeta )}}M(\alpha ,\beta ,\zeta ,t)$$
(77)
we reduce Eq.(75) to the relations which follow from the standard recurrence formulae  for the hypergeometric function $`M(\alpha ,\beta ,\zeta ,t)`$

$`\left(\zeta \alpha \beta +(1t){\displaystyle \frac{d}{dt}}\right)y(\zeta )`$ $`=`$ $`{\displaystyle \frac{(\zeta \alpha )(\zeta \beta )}{\zeta }}y(\zeta +1),`$ (78)
$`\left(\zeta +t{\displaystyle \frac{d}{dt}}\right)y(\zeta +1)`$ $`=`$ $`\zeta y(\zeta ),`$ (79)
where $`y(\zeta )M(\alpha ,\beta ,\zeta ,t)`$.


text:
# Models for the Magnitude-Distribution of Brightest Cluster Galaxies

## 1 Introduction

Among the most luminous bodies in the universe are the brightest, or first-ranked, galaxies in rich clusters. These galaxies have absolute magnitudes between -21.5 and -23.3 and are among the farthest observable objects. In addition, the magnitudes of these brightest cluster galaxies (BCGs) are highly uniform, with a dispersion of 0.32 magnitudes (Hoessel & Schneider 1985). Their uniformity and large luminosity make BCGs excellent standard candles. The uniformity of BCG magnitudes raises a particularly important question regarding their nature (Peebles 1968; Sandage 1972). Are BCGs simply the brightest of a statistical set of galaxies or do they belong to a special class of objects? If a special class of galaxies exists, do all clusters have special galaxies and are they always first-ranked (Bhavsar 1989)? We investigate these questions using extreme value theory (Fisher & Tippett 1928).

## 2 Extreme Value Theory

The motivation for studying extreme phenomena is practical. Many of the memorable experiences in our lives can be classified as statistical extremes. Examples of maximum extremes are floods, the hottest summer temperatures and the lengths of the longest caterpillars. Examples of minimum extremes are draughts, stock market crashes and the wing-spans of the smallest hummingbirds. Some extremes do not effect our lives and others turn them upside down. The desire to understand these types of phenomena prompts the study of extreme value theory.

Fisher & Tippett (1928) show that the distribution of statistically largest or smallest extremes tends asymptotically to a well-determined and analytic form for a general class of parent distributions. Extremes drawn from sufficiently large and steeply falling parent distributions have this form. One may find the original argument in Fisher & Tippett (1928). Their derivation is reconstructed in greater detail by Bhavsar & Barrow (1985), who apply extreme value theory in an analysis of BCG magnitudes. Fisher & Tippett’s result states that the cumulative distribution of maximum extremes is given by:

$$F(x)=e^{e^{a(xx_0)}}.$$
(1)
This distribution is known as the $`Gumbel`$ distribution. (For smallest extremes, one substitutes $`xx`$.) From $`F`$ we may calculate the differential distribution (or probability density):

$$f(x)=ae^{a(xx_0)e^{a(xx_0)}},$$
(2)
where $`f(x)`$ = $`F^{}(x)`$; $`x_0`$ is the mode of the extremes and $`a>0`$ is a measure of the steepness of fall of the parent distribution. The probability density is normalized to unity. The mean, median and standard deviation of the distribution given in Bhavsar & Barrow (1985) correspond to:

$$<x>=x_0+\frac{0.577}{a};med(x)=x_0+\frac{0.367}{a};\sigma ^2=\frac{\pi ^2}{6a^2},$$
(3)
where $`0.577\mathrm{\Gamma }^{}(1)`$ is Euler’s constant, $`0.367ln(ln(2))`$ and $`\sigma `$ is the standard deviation of the extremes. The standard form for the Gumbel, $`F(x)`$ and $`f(x)`$, is shown in Fig. 1. Note that for BCGs, we will be considering minimum extremes (because more negative magnitudes are brighter) and the curves will be inverted ($`xx`$). Henceforth, we will call $`f(x)`$ the Gumbel distribution.

## 3 Brightest Cluster Galaxies

### 3.1 Past Results

Researchers have described BCGs as special, statistical extremes of a normal population and a mixture of the two (Peebles 1968; Peach 1969; Sandage 1972, 1976; Bhavsar & Barrow 1985; Bhavsar 1989; Postman & Lauer 1995).

The motivation for proposing that BCGs are special is due to the small dispersion observed in BCG magnitudes (Peach 1969; Sandage 1972, 1976). These authors argue that such a small dispersion is not sufficiently explained by the steepness of the luminosity function. In addition, astronomers observe a class of BCGs that are morphologically different, called $`cD`$ galaxies. These galaxies are giant ellipticals and often have features, such as multiple nuclei and large envelopes, that distinguish them from normal galaxies.

On the other hand, Peebles (1968) argues that BCGs are just the extreme tail-end of normal galaxies that form in clusters via some stochastic process. In this case, the brightest galaxy in a given cluster is simply the brightest normal galaxy and, therefore, the distribution of BCG magnitudes is a Gumbel. (It is interesting to note that Peebles, independently of Fisher & Tippett, derived the Gumbel distribution for BCGs for the special case of an exponential luminosity function.)

Bhavsar (1989) contends neither of these scenarios adequately describes the observed distribution of BCG magnitudes and argues for a mixed population. Suppose that a special class of Galaxies exists but that not all clusters have a special galaxy. In clusters with no special galaxy, the BCG is simply the brightest normal galaxy. In a cluster containing at least one special galaxy, either all the normal galaxies are fainter, or the brightest normal galaxy(ies) out-shines the special one(s) and attains the first rank. For these reasons, one might expect both types of galaxies to comprise the BCG population. In what follows, we investigate these assumptions quantitatively by analyzing Lauer & Postman’s (1994) data set and revisiting the one used by Bhavsar (1989).

### 3.2 The Distribution of BCG magnitudes

In the case of one population, the distribution function is straight forward. If BCGs are all drawn from a special class of objects, it has been assumed that BCG magnitudes are normally distributed (Peach 1969; Sandage 1972, 1976; Postman & Lauer 1995). In this case, referred to henceforth as model A, the probability distribution of special galaxies, $`f_{sp}`$, is a Gaussian, $`f_g`$, with mean $`M_g`$, standard deviation $`\sigma `$ and normalization such that the integral over all magnitudes, $`M`$, is unity. The distribution function is as follows:

$$f_{sp}(M)=f_g=\frac{1}{\sigma \sqrt{2\pi }}e^{\frac{(MM_g)^2}{2\sigma ^2}}.$$
(4)
If BCGs are simply the brightest of a normal set of galaxies (Peebles 1968), henceforth referred to as model B, the probability distribution of their magnitudes, $`f_{nor}`$, is a Gumbel, $`f_G`$, given by Equation (2), with $`xM`$ and $`x_0M^{}=M_G+\frac{0.577}{a}`$ (Bhavsar & Barrow 1985):

$$f_{nor}(M)=f_G=ae^{a(MM^{})e^{a(MM^{})}},$$
(5)
where $`M_G`$ is the mean of the extremes and $`a`$ is a measure of the steepness of fall of the parent distribution.

In the case of two populations (Bhavsar 1989), we derive the distribution that $`M`$ should have from the contributions of the two individual populations. Consider $`N`$ clusters of galaxies and suppose that $`n<N`$ have at least one special galaxy. Let the independent magnitude-distribution of normal and special galaxies, respectively, be $`f_{nor}`$ and $`f_{sp}`$. The total magnitude-distribution function, $`f_{tot}`$, is then given by:

$`f_{tot}(M)`$ $`=`$ $`d[f_{sp}{\displaystyle _M^{\mathrm{}}}f_{nor}(M^{})dM^{}+`$ (6)
$`f_{nor}{\displaystyle _M^{\mathrm{}}}f_{sp}(M^{})dM^{}]+`$
$`(1d)f_{nor},`$
where $`d=n/N`$. The first (second) term is the probability of picking a special (normal) galaxy, with absolute magnitude $`M`$, from a cluster containing both populations with the condition that all the normal (special) galaxies are fainter. The third term gives the probability of picking a galaxy, with absolute magnitude $`M`$, in clusters containing only normal galaxies. Equation (6) is true for all well-behaved functions $`f_{nor}`$ and $`f_{sp}`$. If $`f_{nor}`$ and $`f_{sp}`$ are normalized to unity, then so is the resulting total distribution function $`f_{tot}`$. (Note that Equation (6) works, in general, whenever there are two independent populations competing for first rank.)

For BCGs, we consider three different two-population models. The first is the case discussed above with the brightest normal galaxies comprising one population and a special class of galaxies comprising the other. We call this case model C and write the total distribution as $`f_{Gg}`$ (where ‘$`Gg`$’ stands for ‘$`G`$umbel + $`g`$aussian’). To obtain the final form of $`f_{Gg}`$ we note that:

$$I_G=_M^{\mathrm{}}f_G(M^{})𝑑M^{}=F(M),$$
(7)
where $`F(M)`$ is given by Equation (1) with $`xM`$, and $`x_0M_G+\frac{0.577}{a}`$. Second, we note that:

$$I_g=_M^{\mathrm{}}f_g(M^{})𝑑M^{}=(1\pm erf|MM_g|)/2,$$
(8)
where $`erf`$ is the error function. The upper sign is for $`M<M_g`$ and the lower sign is for $`M>M_g`$. Thus, we may rewrite $`f_{Gg}`$ by substituting in $`I_G`$ and $`I_g`$:

$$f_{Gg}(M)=d[f_gI_G+f_GI_g]+(1d)f_G.$$
(9)
Other possible combinations of assigning $`f_G`$ and $`f_g`$ to the two populations result in models D and E. In the case of model D, both distributions are Gaussian and the total distribution function, $`f_{gg}`$, is given by:

$$f_{gg}(M)=d[f_{g2}I_{g1}+f_{g1}I_{g2}]+(1d)f_{g1},$$
(10)
where the notation is self-evident and the two Gaussians are characterized, respectively, by $`M_{g1}`$, $`\sigma _1`$ and $`M_{g2}`$, $`\sigma _2`$. In the case of model E, both distributions are Gumbels ($`f_{sp}`$ is also a Gumbel) and the total distribution function, $`f_{GG}`$, is given by:

$$f_{GG}(M)=d[f_{G2}I_{G1}+f_{G1}I_{G2}]+(1d)f_{G1},$$
(11)
where the two Gumbels are characterized, respectively, by $`M_{G1}`$, $`a_1`$ and $`M_{G2}`$, $`a_2`$. Table 1 summarizes the forms of the five models.

## 4 Modeling the Data

### 4.1 Data Sets

We utilize two data sets from the literature. First, we reanalyze the data used by Bhavsar (1989). This is a 93 member subset of 116 metric BCG visual-intrinsic (VI) magnitudes compiled by Hoessel, Gunn & Thuan (1980), henceforth referred to as “HGT”. These 93 are the data from clusters of richness 0 and 1 only; Bhavsar ignores the rest of the BCGs in order to keep the data set homogeneous. The BCG magnitudes are internally consistent to 0.04 magnitudes, as published in HGT. Second, we analyze the 119 metric BCG magnitudes, taken in the Kron-Cousins $`R_c`$ band, compiled by Lauer & Postman (1994), henceforth referred to as “LP”. The data were corrected for local and possible large scale galactic motions. The 119 LP data are comprised of BGCs from 107 clusters of richness 0 & 1, and 9, 2 and 1 of richness 2, 3 and 4, respectively (Abell, Corwin & Olowin 1989). We find that removing the 12 BCGs from clusters of richness class $``$ 2 does not significantly change the distribution of the LP data. This is consistent with Sandage’s (1976) result that BCG magnitude is independent of cluster-richness. The internal consistency of the set is 0.014 magnitudes, as published in Postman & Lauer (1995). Bhavsar (1989) proposes a two-population model for the HGT data. His maximum-likelihood fit is consistent with the data and has parameter-values consistent with physically measured quantities. Postman & Lauer (1995) conclude that the LP data are consistent with a Gaussian.

There are differences in the data sets that could be the reason for the disagreement between Bhavsar (1989) and Postman & Lauer (1995). The two were obtained in different optical bands. The mean of the HGT data set is 0.2 magnitudes brighter than the mean of the LP data set. The two data sets have 34 galaxies in common. Comparing the subset of 34, we find that the HGT values are, on average, 0.06 $`\pm `$ 0.19 magnitudes brighter than the LP values. A two-sample Kolmogorov-Smirnov (K-S) Test addresses the consistency of the two data sets in describing the same population of objects. The null hypothesis is that the same distribution describes both data sets. We find that the two data sets fail the null hypothesis at the 82% confidence level. Therefore, we do not expect the same parameters or distribution to describe both sets. These discrepancies may need further investigation, but such an analysis is outside the scope of this work. We investigate each data set separately and present our results.

### 4.2 Fitting Method

We consider models A-E discussed above. The two-population distributions have five parameters each: two means, two standard deviations and the fraction, $`d`$, of clusters that contain a special population of galaxies. If there is no population of special galaxies, then $`d`$ = 0. We use maximum-likelihood fitting. The theory behind this method is discussed in Press, et al. (1992). The Maximum-Likelihood fit to a data set of size $`N`$ for a function, $`f`$, are the parameters, $`𝐚`$, that maximize the likelihood function:

$$L=\underset{i=1}{\overset{N}{}}f(x_i;𝐚),$$
(12)
where the $`f(x_i;𝐚)`$ are the values of the probability density, $`f`$, evaluated at each of the $`N`$ data points, $`x_i`$. For a certain $`f`$, one finds the set of parameters that maximizes the product, $`L`$.

## 5 Results

### 5.1 Parameters and Fits

After obtaining parameters from the maximum-likelihood method for models A-E for both data sets, we compute the K-S statistics. We list the results in Tables 2 & 3, respectively. Lower values of the K-S $`D`$-statistic correspond to lower values of rejection probability, $`P`$, and thus denote a better fit. Figs. 2 & 3 illustrate the performance of each of the five models. Note that the distributions use the parameters obtained by the maximum-likelihood method, using $`every`$ data point, and are $`not`$ a fit to the particular histograms.

### 5.2 Comparison with Previous Work

We compare our results with Bhavsar (1989) and Postman & Lauer (1995). Bhavsar’s (1989) two-population model is our model C. He uses maximum-likelihood fitting and his best-fitting parameters are $`M_G`$ = -22.31, $`M_g`$ = -22.79, $`d`$ = 0.63, $`a`$ = 4.01 and $`\sigma `$ = 0.21. Our parameters are in excellent agreement. Minor variation is expected due to differences in fitting techniques. Postman & Lauer (1995) argue against Bhavsar’s two-population model and claim that BCG magnitudes are Gaussian, based on a 26% confidence level.

In agreement with both Bhavsar (1989) and Postman & Lauer (1995), it is clear from Tables 2 & 3 and Figs. 2 & 3 that for both data sets no Gumbel distribution describes the BCG data. This rejects the Gumbel hypothesis (model B) with 85% and 93% confidence levels, respectively, for the HGT and LP sets. For the HGT data, the Gaussian fails at the 53% confidence level, while for the LP data, the rejection confidence is 16%. The difference between our value of 16% and Postman & Lauer’s value arises because our result is for the maximum-likelihood fit Gaussian, while Postman & Lauer’s is for a Gaussian with the same mean and standard deviation as the LP data.

The relatively high rejection-confidence of the one-population models has motivated us to investigate two-population models. The presence of cD galaxies strongly suggests the possibility of another population. Overall, the two-population models fit the data much better than do the Gumbels and as well or better than do the respective Gaussians. The larger number of parameters is taken into account by the statistical estimators when calculating the confidence of rejecting the null hypothesis. Moreover, the parameters are physical quantities that are observationally verifiable (Bhavsar 1989).

Our result that no one model or set of parameters describes both data sets is consistent with the fact that a two-sample K-S Test indicates that the sets are not consistent with one another. Postman & Lauer (1995) have raised questions regarding HGT’s BCG classification and sky subtraction.

### 5.3 Physical Motivation

Researchers have suggested various mechanisms whereby a second population with a brighter average metric magnitude could evolve from the bright normal galaxies. Cannibalism, the process by which large galaxies in the central regions of rich clusters grow at the expense of smaller galaxies (Ostriker & Hausman 1977; Hausman & Ostriker 1978), is one possibility. The existence of giant elliptical and cD galaxies near the centre of approximately half of all rich clusters supports this hypothesis. These galaxies always lie at the tail-end of their cluster-luminosity functions. The occurrence of cannibalism continues to be debated (Merritt 1984).

Motivated by the existence in the literature of strong arguments for such a process, we build a very simple schematic to study its $`statistical`$ effects on the population of first-ranked galaxies. We make two assumptions: (i) at an early epoch the BCGs all belonged to one population and (ii) galaxies from the bright end of this population evolve, resulting in a random boost to their luminosity. We construct a set of $`N`$ galaxies with an exponential luminosity function between absolute magnitudes -22.0 and -23.0. This represents the galaxies at the bright end of cluster luminosity functions that are candidates for a boost. A random number, $`n`$, of these galaxies undergo a random boost between 0.1 and 0.9 magnitudes. We label the boosted subset as $`n_b`$. We choose this range for the following reasons. First, Hausman & Ostriker (1978) show via a simulation that one would expect a large galaxy to gain, on average, 0.5 magnitudes during its first cannibalistic encounter. This is consistent with Aragon-Salamanca, Baugh & Kauffmann (1998), who state that BCGs were approximately 0.5 magnitudes fainter at $`z`$ = 1. Second, we limit ourselves to one encounter because Merritt (1984) argues that the time scale for galactic encounters is too long for cannibalism to be common in the universe. We wish to investigate the magnitude-distribution of the resulting boosted population. These represent the special galaxies mentioned previously. Specifically, this distribution could give us insight into the form of $`f_{sp}`$.

To our surprise, we find that the distribution, $`f_{sp}`$, of $`n_b`$ is a Gumbel! The K-S Test rejects the Gaussian hypothesis at the 98% confidence level. Conversely, the Gumbel distribution, with the same mean and deviation as the data, fits well, with only a 7% confidence level for rejection. We summarize these results in Table and Fig. 4. Thus, the two-population model E (a combination of two Gumbels), which is best-fitting for the newer LP data, has a physical basis.

## 6 Conclusion

For more than thirty years, cosmologists have debated the nature of the magnitude-distribution of brightest cluster galaxies. Peebles (1968) and Sandage (1972, 1976) & Peach (1969) reach markedly different conclusions. More recently, Bhavsar (1989) and Postman & Lauer (1995) differ regarding the population(s) that comprise the first-ranked galaxies. In light of this controversy, we have conducted a new examination of the distribution of BCG magnitudes. We consider the BCGs as $`a`$ $`class`$ $`of`$ $`objects`$ to which we may apply well established results from extreme value theory. We find that there are a number of models that perform well in describing the HGT and LP data sets. Though a Gaussian fits both data sets, the confidence limits warrant further investigation of two-population models.

Tables 2 & 3 clearly show that we should reject the Gumbel (model B) as a fit, i.e., the hypothesis that all BCGs are statistical extremes. The Gaussian (model A) is marginally acceptable but without physical basis. Two-population models, in particular, the three combinations of $`f_G`$ and $`f_g`$, describe the data very well. Tables 2 & 3 show their relative merits. Model E stands out as giving the best overall fit and is motivated by a physical basis. Therefore, it is most likely that there are two populations of BCGs: the extremes of a normal population and a class of atypical galaxies with a brighter average mean.

We thank Marc Postman for sending us the LP data. This research was supported by an ANN grant from the US Department of Education and the Kentucky Space Grant Consortium.