text:
# On thermostats and entropy production

## I Introduction

There has been an intense interest over the past several years in the properties of fluid systems subjected to Gaussian thermostats***In fact Nosé-Hoover thermostats have been the preferred tool for doing simulations at a given temperature. These share all properties discussed here for Gaussian thermostats.. These thermostats are fictitious reversible force fields added to the equations of motion for fluid systems, and were introduced by people doing molecular dynamics simulations on fluids in stationary nonequilibrium systems in order to prevent the fluids from heating up indefinitely. These thermostatted systems have attracted the attention of workers in statistical mechanics and in dynamical systems theory for a variety of reasons: Firstly, since phase space volumes are not conserved in systems with Gaussian thermostats, and since the irreversible entropy production in systems subjected to these type of thermostats can be equated to the average rate of contraction of phase space volume on a microscopic scale, it follows that the rate of entropy production is equal to the negative of the sum of all Lyapunov exponents in the system . Secondly, if, as often is the case, the irreversible entropy production may also be expressed in terms of transport coefficients and thermodynamic driving forces, there is a direct link between dynamical systems properties and nonequilibrium statistical mechanics for such systems. Thirdly, a point of great interest is that the phase space contraction induced by Gaussian thermostats in macroscopically stationary nonequilibrium states causes these states to live on fractal attractors that can be characterized by Sinai-Ruelle-Bowen (SRB) measures. Gallavotti, Cohen and Ruelle have conjectured that these may be the generic type of states needed to describe stationary nonequilibrium states. Their Kaplan-Yorke dimension, which is related directly to sums of Lyapunov exponents, also measures the irreversible entropy production as recently discussed by Evans et  al., in a way similar to that described by Gaspard for Hamiltonian systems with escape. Finally, Ruelle has argued that in a stationary state on a compact subset of phase space, such as a shell of constant kinetic energy, the average rate of phase space contraction cannot be negative, since the phase space volume occupied by the state cannot expand without bounds. Therefore the irreversible entropy production cannot be negative, which would demonstrate the validity of the Second Law of thermodynamics for Gaussian thermostated systems.

In order to avoid having fictitious Gaussian thermostats acting everywhere in the system, Chernov and Lebowitz studied a different class of thermostats that are active only at the boundaries of the system. These thermostats act in such a way that no energy is exchanged between thermostat and system, but the velocity distribution of the particles reflected from the walls of the system is contracted in such a way that a stationary shear flow results, giving rise to an irreversible entropy production. In the infinite system limit the shear induced irreversible entropy production could be equated again to the average rate of phase space contraction produced by the thermostat, but other cases no such relation existed for finite systems clear deviations from this equality were found. Obviously this could be due either to irreversible entropy production in kinetic modes, which are always present in boundary layers, or to a a complete breakdown of the relation between phase space contraction rate and irreversible entropy production. Related results have also been obtained by Klages and coworkers. These authors have constructed a variety of thermostats which act only during collisions of particles with fixed scatterers or walls, and also lead to nonequilibrium steady states, with phase space contraction onto a lower dimensional attractor. In these models, the entropy production may or may not be related to the rate of phase space contraction, depending upon the precise nature of the thermostat and the system.

Here we will show on the basis of a simple and fairly natural example that indeed one can find thermostats for which the relation between the rate of phase space contraction and the rate of entropy production breaks down, and, in addition, we show one may have stationary nonequilibrium states without an overall phase space contraction, as previously shown by Gaspard. Our model consists of a Lorentz gas enclosed between walls at different temperatures, a system for which Spohn and Lebowitz proved long ago that it satisfies Fourier’s law of heat conduction in the Grad limit, in which one sends the size of the scatterers to zero while increasing their density in such a way that the ratio between mean free path and system size remains constant. The support of the stationary measure will have the same dimension as the full phase space, yet it will exhibit strong fractal features in the form of a fractal repeller located on the boundary between two subsets of different energy, very similar to the fractal boundaries studied by Gaspard in the framework of particle diffusion . The Lyapunov exponents of the system will turn out to be simply related to those of the equilibrium system. Finally we will comment on the consequences of our considerations for the Second Law. We will argue that the assumptions used to obtain this law in our context may be replaced by other, equally reasonable, assumptions which lead to violations of the Second Law.

## II The Lorentz gas

The system considered here is the usual Lorentz gas consisting of fixed spherical scatterers located at random positions in a $`d`$-dimensional space, and a light point particle moving between the scatterers at constant speed and making specular, elastic collisions with them. We will assume that the scatterers cannot overlap each other so as to avoid some problems with ergodicity and percolation that may occur in the case of overlapping scatterers. Instead of this choice we could also consider periodic billiards, with scatterers arranged on a regular lattice (preferably without infinite horizons), or systems with randomly oriented polyhedral scatterers instead of spherical ones. In the latter case all Lyapunov exponents would be trivially zero, but, as in the case of the periodic billiards, the transport properties would hardly differ from those of the ordinary Lorentz gas. As boundaries we choose two parallel flat walls, separated by a distance $`L`$. The system may be infinite or periodically repeated in the directions parallel to the walls. The collisions of the light particle with these walls are also specular in direction, but the speed with which the particle leaves a wall will always be uniformly distributed in the range $`v_i(1ϵ/2)vv_i(1+ϵ/2),i=1,2,`$ where $`ϵ<<1`$, and $`v_1=(dk_BT_1/m)^{1/2}`$ is the average speed for one wall and $`v_2=(dk_BT_2/m)^{1/2}`$, for the otherTo be slightly more realistic one could sample $`v_1`$ and $`v_2`$ from Maxwellian distributions with temperatures $`T_1`$ respectively $`T_2`$, but for our arguments this would not make any difference.. Eventually we have in mind the limiting case where $`ϵ`$ tends to zero. Further, $`m`$ is the mass of the light particle, $`T`$ is temperature, $`k_B`$ Boltzmann’s constant and $`d`$ the dimensionality. We suppose that the scatterers are distributed in space with number density $`n`$, and that they have radius $`a`$.

Spohn and Lebowitz have shown rigorously that in the Grad limit, $`n\mathrm{},a0,`$ with $`0<na^{d1}<<1`$, combined with the limit of $`L`$ tending to infinity, this system satisfies Fourier’s law of heat conduction: the temperature profile becomes linear and the energy current through the system becomes proportional to the temperature gradient. In addition, the coefficient of heat conduction can be obtained from the Lorentz-Boltzmann equation in this limit. One expects these results to hold much more generally in the limit of large L, except that for higher densities of scatterers (measured in dimensionless units $`\stackrel{~}{n}=na^d`$, with $`\stackrel{~}{n}>0`$) the transport coefficient will not be following from the Lorentz-Boltzmann equation any more.

## III Entropy production and dynamical properties

It is easy to find an expression for the irreversible entropy production for these systems. As is usual for stationary thermostated systems one has to argue that the irreversible entropy production per unit time has to equal the rate of entropy flow from the system into the thermostats. In the stationary state the average energy current through the system has to have a constant value $`j_e`$, so the average total entropy production per unit time $`\sigma `$ satisfies

$$\sigma =j_e(\frac{1}{T_2}\frac{1}{T_1}),$$
(1)
where we assume $`T_1>T_2`$. In the present model it follows immediately that the energy current is directed from higher to lower temperature, as every time a particle changes its velocity from $`v_1`$ to $`v_2`$ it transfers an amount of energy $`dk_B(T_1T_2)/2`$ to the cold reservoir, whereas in changes from $`v_2`$ to $`v_1`$ it extracts the same amount of energy from the hot reservoir. It is not hard making explicit estimates of $`j_e`$, especially for small values of $`\stackrel{~}{n}`$, but for our present purposes it suffices that this current is non-zero. The phase space density of the particles traveling in the small range of speeds about $`v_i`$ is easily seen to be proportional to $`v_i^{(d+1)}`$. This ensures, among other things, that the number of particles incident on a small area of the wall per unit time is equal to the number of particles leaving the small area per unit time. Furthermore, under the action of the thermostats as defined above there is no net contraction of phase space density. Indeed on a collision reducing $`v_1`$ to $`v_2`$ the phase space density at the coordinates of the collision is multiplied by a factor $`(v_1/v_2)^{(d+1)}`$, but on the next collision bringing the speed back to $`v_1`$ this contraction is undone exactly and therefore the average rate of phase space contraction over a long time becomes exactly zero. As a result the identity between average phase space contraction rate and irreversible entropy production that holds for Gaussian thermostats, is not valid in the present case. In addition, the presence of a stationary heat current influences the Lyapunov exponents in a trivial way only. The crucial observation here is to note that the trajectories of the light particle are exactly the same as in equilibrium, only the parts of it that originate from the first wall are traversed at speed $`v_1`$ and those originating from the second wall with speed $`v_2`$. As a result the total time required for traversing a long trajectory is multiplied by a factor of $`(1/v_1+1/v_2)/(2/v)`$ compared to the time needed by a particle with constant speed $`v`$. Hence all the Lyapunov exponents are multiplied by the inverse of this factor. The stationary measure on phase space is a measure with one constant density, proportional to $`(1/v_1)^{(d+1)}`$ and restricted to the velocity range $`v_1(1ϵ/2)vv_1(1+ϵ/2)`$, for points that can be traced back to a last collision with the first wall, and another constant density, proportional to $`(1/v_2)^{(d+1)}`$ and confined to the velocity range $`v_2(1ϵ/2)vv_2(1+ϵ/2)`$, for points coming from a last collision with the other wall. It is completely analogous to the measures considered by Gaspard to describe systems with a stationary diffusion current.

In spite of the absence of phase space contraction and the simplicity of the Lyapunov spectrum these stationary measures have strong fractal properties. Due to the way the model is constructed, the phase space trajectories for the moving particles are located on two thin energy shells about energies $`mv_i^2/2,i=1,2`$. In the stationary state, the trajectories on each shell are those whose last collision with a wall led to a value of the energy in the range of the shell. The trajectory of the particle may, of course, take place on both shells. To better visualize this phase space structure, it is helpful to rescale each of the energy shells so that they coincide. Under these circumstances, one may visualize the stationary state phase space regions as consisting of two disjoint subsets, separated by a boundary. Each of the subsets consists of points whose motion, extrapolated back in time reaches a particular wall. The boundary between these two subsets consists of points in phase space from which, extrapolating back in time, neither boundary is ever reached. An even smaller subset of this consists of points from which the boundaries are never reached either in the forward or in the backward time direction. This latter set is called a repeller, and it coincides with the intersection of its stable and unstable manifolds, that is, with the intersection of the set of points which may have collided with one or the other walls in the past but will not do so in the future (the stable manifold of the repeller), or have not collided with a wall in the past but may do so in the future (the unstable manifold of the repeller). Therefore, it is clear that the boundary of the two sets under discussion here coincides with the unstable manifold of the repeller. The fractal repeller may contain subsets of trajectories with an environment extrapolating back entirely to a single boundary. These will not appear in the boundary set considered presently. However, we expect that typical trajectories on the repeller will extend throughout the system and get arbitrarily close to either boundary. Hence the missing part should be of vanishing measure relative to the main body of the repeller.

The full fractal repeller described above plays a central role in the escape rate formalism of Gaspard and Nicolis. In this formalism the escape rate from a system with open boundaries is related on the one hand to the diffusion coefficient, $`D`$, and on the other hand to the sum of positive Lyapunov exponents, $`\lambda _i^+`$ and the Kolmogorov-Sinai (KS) entropy, $`h_{KS}`$, on the repeller. For a Lorentz gas in two dimensions, Gaspard has shown how to express the partial information dimension, $`d_u`$, of the unstable manifold of the repeller to the diffusion coefficient of the moving particle in the fixed array of scatterers, and to the Lyapunov exponents and KS entropy. In our case, this fractal dimension is given by

$`d_u`$ $`=`$ $`{\displaystyle \frac{h_{KS}}{\lambda ^+}},`$ (2)
$`=`$ $`1{\displaystyle \frac{D\pi ^2}{\lambda ^+L^2}}+o(L^2).`$ (3)
Simple scaling arguments show that this dimension is independent of the velocity of the particle on the unstable manifold, so that our scaling of the two energy shells to make them coincide does not affect the dimension of the boundary separating the two sets of phase points.

In the present system the coefficient of heat conduction may also be related to properties of the repeller. This follows from the work of Dorfman and Gaspard who derived a relation between the coefficient of heat conductivity and the dynamical properties of trajectories on an appropriate repeller for heat conduction. This repeller is obtained by looking at the fluctuations of the Helfand moments associated with heat conduction. The Helfand moments undergo a deterministic diffusion in phase space. When the magnitude of the Helfand moment for a given trajectory in phase space exceeds a certain value, say $`(d/2)k_BT|L/2|`$, one considers that the system has “escaped” from a region in phase space, and the escape-rate formalism may be applied. In the case of the Lorentz gas under discussion here, the Helfand moments for heat flow and for diffusion are essentially the same, differing only by a constant factor, since the moving particle keeps a constant energy one collision with a wall to the next. until collIding with a wall. Thus the dimensional properties of the stable and unstable manifolds of the repellers are the same in both cases, diffusion and heat conductivity.

## IV Discussion

The simple example discussed in this paper shows that macroscopically stationary nonequilibrium measures need not be phase space contracting and, consequentially, the average rate of irreversible entropy production need not be equal to the average rate of phase space contraction. In fact it seems to us that the thermostat considered here is closer to physical reality than the Gaussian thermostats that do give rise to the identities mentioned above. The net effect of our thermostat may be described as mimicking a transport of phase space density between the hot and the cold reservoir. As a result of this the temperatures of these reservoirs will approach each other (though very slowly if the reservoirs are really large) and the coarse grained entropy, based on a coarse grained phase space density, will decrease in the hot reservoir, but increase by a larger amount in the cold reservoir. This scenario can be made very plausible by probabilistic arguments, but is very hard to prove for almost any Hamiltonian model. Both Gaussian thermostats and the thermostat considered above can be considered as simplified models of real, Hamiltonian thermostats, in which the expected properties have been built in from scratch. It does not appear plausible though that in realistic systems coupled to thermostats at different temperatures the microscopic phase space density will keep increasing forever without bounds. In view of Liouville’s theorem physical changes in phase space density always must be the result of exchanges of phase space density between the system and its thermostats. If the latter may be treated as stationary sources the distribution of phase space density in a stationary nonequilibrium state for the system has to be stationary as well. In our opinion the feature of steady phase space contraction really is an artifact of Gaussian thermostats. This even raises the question whether there is a real use for these besides their handiness in nonequilibrium steady state simulations, but there certainly is. First of all Gaussian thermostats can be used as dynamical probes of nonequilibrium properties, precisely because of their satisfying the equality of entropy production and phase space contraction rate. Further, there are strong indications that at least for systems not too far from equilibrium, Gaussian thermostats lead to behavior that is almost indistinguishable from that brought about by more physical thermostats. The nonequilibrium features of stationary nonequilibrium measures show up predominantly in their fractal features. The latter are much more visible for systems with Gaussian thermostats, in which the stationary measure lives on a fractal exclusively, than in the presumably more physical stationary measures in which the fractal features only show up in the boundaries between subsets with different densities. Therefore systems with Gaussian thermostats may be more suitable in the end for a dynamical study of nonequilibrium states. At this moment it is unclear whether this really is the case and in addition the relations between systems with different types of thermostats still largely have to be established.

We now make some remarks on the implications of thermostating mechanisms for the Second Law of thermodynamics. We indicated already that the generation of an average energy current from the hot to the cold reservoir has been built into our thermostat by hand. To make this point somewhat clearer, let us consider a somewhat more physical thermostat. As illustrated in Figure 1 it consists of a gas of particles that do not see the Lorentz gas walls, but instead are confined by a wall somewhat inside the Lorentz gas region. The latter wall in turn is invisible to the light Lorentz gas particle. In the narrow strip between these two walls the Lorentz gas particle may collide with the thermostat particles and so, if these collisions are sufficiently frequent, the light particle will leave the strip with a kinetic energy that on average is proportional to the temperature of the thermostat.

The combined action of two of these thermostats on either side, with temperatures $`T_1`$ and $`T_2`$, will be basically the same as that of our simple model thermostats. For this to be true, however, we implicitly have to break time reversal symmetry by assuming that the velocity distributions of the Lorentz particle and of the thermostat particles before their collision are mostly uncorrelated, in other words we have to assume Boltzmann’s Stoszzahlansatz or some generalization of it. Instead, we could have required such property to hold after the collision . This would correspond to the macroscopic arrow of time running in the opposite direction. Of course, we as observers would not notice the difference. Only an outside observer could notice our inversion of the notions of past and future. But, as remarked by e.g. Schulman, more exotic possibilities cannot be excluded. The arrow of time could run one way in parts of the system and the other way in other parts. Also, on the basis of known microscopic equations of motion one cannot exclude the possibility that at some time the direction of the macroscopic arrow of time will turn around, perhaps through growth of the area in which it runs the opposite way. To explain our obvious observation of a well defined direction of time, reversible microscopic equations simply do not suffice. Apparently we do need the additional fact that at some point in time our universe was brought into a state of extremely low entropy, so we observe an arrow of time moving away quite definitely from that point. And in spite of any lack of evidence we cannot exclude the possibility that the temporal boundary conditions to the ‘equations of motion of our universe’ are not all located in the past, which then might lead indeed to eventual inversions of the arrow of time. In view of this it seems clear that proofs of the Second Law will have to require additional information about the properties and the history of our universe and cannot be based exclusively upon the mixing or chaotic properties of the microscopic equations of motion, or on properties of thermostats in stationary nonequilibrium systems. In less cosmic language, we might say that some assumptions about the absence of initial correlations in our systems will have to be made in order to demonstrate the validity of the Second Law, no matter how chaotic the dynamics of our system. Along these lines, a much more modest goal, which is still very hard to attain, is to show that isolated systems showing sufficient lack of correlations initially, will indeed satisfy the Second Law. Here perhaps thermostats can be of some use, but in order to be credible, they would have to be of Hamiltonian type .

## V Acknowledgements

For both of us Joel Lebowitz has been a dear friend for over thirty five years. It is a great pleasure being able to contribute to this very special Festschrift. We congratulate Joel on his seventieth birthday and express our hope he will continue for many more years being an inspiring example for all statistical physicists in all his scientific and humanitarian efforts. We thank Rainer Klages for some valuable remarks on the manuscript. JRD wishes to thank FOM for support of visits to the University of Utrecht, and the National Science Foundation (USA) for support under grants PHY-96-00428 and PHY-98-20824. HvB acknowledges support by FOM, SMC and by the NWO Priority Program Non-Linear Systems, which are financially supported by the ”Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)”.


format: text

text:
# Stability criterion for multi-component solitary waves

## I Introduction

Recent discovery of self-focusing of partially coherent light and experimental observation of the so-called incoherent spatial solitons  has called for a systematic analysis of the properties and stability of multi-component and multi-parameter solitary waves. Incoherent solitons are generated in noninstantaneous nonlinear media such as biased photorefractive crystals. In this case, a self-consistent modal theory , which is equivalent to the coherent density approach, describes the incoherent solitons with the help of a system of coupled nonlinear Schrödinger (NLS) equations (see also ). Similar models appear, in a different physical context, in the theory of soliton wavelength-division multiplexing , multi-channel bit-parallel-wavelength optical fiber networks , multi-species and spinor Bose-Einstein condensates , and other important applications . In all such physical models solitary waves are multi-component, being described by localized solutions of the coupled nonlinear equations. In some very special cases, the coupled system allows for explicit analytical solutions (see, e.g., Ref. ) but, generally speaking, the nonlinear models with multi-component solitary waves are nonintegrable. Stability of solitary waves is therefore a crucial issue for any kind of their applications.

The study of soliton stability has a long history. Stability of one-parameter solitary waves has been already well understood for both fundamental (single-hump and nodeless) solitons  and solitons with nodes and multiple humps . The pioneering results of Vakhitov and Kolokolov  found their rigorous justification in a general mathematical theory of Grillakis, Shatah, and Strauss . Although the corresponding stability and instability theorems for the scalar NLS models extend formally to the case of multi-parametric solitons , all the cases analyzed so far correspond to solitary waves with effectively a single parameter.

A recent progress in the study of soliton instabilities is associated with the application of a bifurcation theory valid for weakly unstable stationary localized waves. In this method, the corresponding unstable eigenvalue of the associated spectral problem is treated as a small parameter of multi-scale asymptotic expansions . In the case of multi-parameter solitary waves, a simplified version of this method is usually reduced to a number of “magic determinants” constructed from the derivatives of the system invariants near a marginal stability line . However, such a bifurcation method has no rigorous proof, and it does not allow to predict the complete domains of the soliton stability and instability, since more general oscillatory instabilities may occur as well .

In this paper, we present a complete theory for stability and instability of multi-parameter solitary waves considering a particular example of a system of $`N`$ incoherently coupled NLS equations. Our results include the asymptotic bifurcation method with the determinant criterion as a simple near-threshold limiting case. They also expand the applicability boundaries of the previously known mathematical theorems  to the case of multi-component and multi-parameter solitary waves.

The system of incoherently coupled NLS equations has already been studied in many papers (see, e.g., Refs.  to cite a few). However, the study of stability of single-hump and multi-hump solitary waves was restricted again by a single-parameter case, when soliton components have a similar shape and their amplitudes are proportional to each other . In this paper, we expand those results and present, for the first time to our knowledge, a complete matrix analysis of the constrained variational problem leading to the finite-dimensional linear algebra. Although some of our results depend on the properties which are specific to the model under consideration, we believe that both the method and the basic results can be generalized, under proper assumptions, to be applied to other types of nonlinear physical models that support multi-parameter solitary waves.

## II Model and Basic Results

We consider nonlinear interaction of $`N`$ optical modes that describe either a multi-mode structure of a partially incoherent self-trapped beam or incoherent coupling between optical channels with different wavelengths in a fiber. Then, the amplitude envelopes of partial modes satisfy the following system of incoherently coupled NLS equations,

$$i\frac{\psi _n}{z}+d_n_𝐱^2\psi _n+\left(\underset{m=1}{\overset{N}{}}\gamma _{nm}|\psi _m|^2\right)\psi _n=0,$$
(1)
where $`_𝐱^2`$ stands for Laplacian in the $`D`$dimensional space $`𝐱=(x_1,\mathrm{},x_D)`$, and all the coefficients $`d_n`$ are assumed to be positive. When one of the variables of the vector $`𝐱`$ stands for time, Eqs. (1) describe the spatio-temporal dynamics of self-focused and self-modulated light in the form of the so-called light bullets.

Provided the symmetry conditions $`\gamma _{nm}=\gamma _{mn}`$ are satisfied, the system (1) conserves the Hamiltonian,

$`H={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑𝐱\left({\displaystyle \underset{n=1}{\overset{N}{}}}d_n\left|_𝐱\psi _n\right|^2{\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{N}{}}}{\displaystyle \underset{m=1}{\overset{N}{}}}\gamma _{nm}|\psi _n|^2|\psi _m|^2\right),`$

the individual mode powers, $`Q_n=\frac{1}{2}|\psi _n|^2𝑑𝐱`$, and the total field momentum. Localized solutions of Eqs. (1) for the fundamental solitary waves are defined as $`\psi _n=\mathrm{\Phi }_n(𝐱)e^{i\beta _nz}`$, where $`\mathrm{\Phi }_n(𝐱)`$ are real functions with no nodes, and $`\beta _n`$ are positive propagation constants. The soliton solutions are stationary points of the Lyapunov functional,

$$\mathrm{\Lambda }[𝝍]=H[𝝍]+\underset{n=1}{\overset{N}{}}\beta _nQ_n[𝝍],$$
(2)
i.e. the first variation of $`\mathrm{\Lambda }[𝝍]`$ vanishes at $`𝝍=𝚽(𝐱)`$. The second variation of $`\mathrm{\Lambda }[𝝍]`$ defines the stability properties: negative directions of the second variation corresponds to unstable eigenvalues in the soliton stability problem (see, e.g., Ref.  for a review of the basic results).

The stability problem is defined by minimizing the second variation of the Lyapunov functional $`\mathrm{\Lambda }[𝝍]`$,

$$\delta ^2\mathrm{\Lambda }=_{\mathrm{}}^{\mathrm{}}𝑑𝐱\left[𝐮|𝐋_1𝐮+𝐰|𝐋_0𝐰\right],$$
(3)
where $`𝐮(𝐱)`$ and $`𝐰(𝐱)`$ are perturbations of the multi-component solitary wave taken in the form $`𝝍=𝚽(𝐱)+\left[𝐮+i𝐰\right](𝐱)e^{\lambda z}`$, and the scalar product is defined as $`𝐟|𝐠=_{n=1}^Nf_n^{}g_n`$. The matrix Sturm-Liouville operator $`𝐋_0`$ has a diagonal form with the elements

$$(L_0)_{nn}=d_n_𝐱^2+\beta _n\underset{m=1}{\overset{N}{}}\gamma _{nm}\mathrm{\Phi }_m^2,$$

and the matrix operator $`𝐋_1`$ has the elements

$$(L_1)_{nn}=d_n_𝐱^2+\beta _n\underset{m=1}{\overset{N}{}}\gamma _{nm}\mathrm{\Phi }_m^22\gamma _{nn}\mathrm{\Phi }_n^2,$$

at the diagonal, and $`(L_1)_{nm}=2\gamma _{nm}\mathrm{\Phi }_n\mathrm{\Phi }_m`$, off the diagonal. Operators $`L_0`$ and $`L_1`$ determine the linear eigenvalue problem for stability of multi-component solitary waves,

$$𝐋_1𝐮=\lambda 𝐰,𝐋_0𝐰=\lambda 𝐮.$$
(4)
Both the linear problem (4) and minimization problem (3) should satisfy a set of $`N`$ constraints,

$$F_n=_{\mathrm{}}^{\mathrm{}}𝑑𝐱\mathrm{\Phi }_n𝐞_n|𝐮=0,$$
(5)
where $`𝐞_n`$ is the $`n^{\text{th}}`$ unit vector, which correspond to the conservation of the individual powers $`Q_n`$ under the action of a perturbation described by a vector ($`𝐮,𝐰`$).

First of all, we recall the main result of Refs.  that the one-parameter solitary waves with no nodes ($`N=1`$) are stable in the framework of the constrained variational problem (3)-(5) provided the energetic surface $`\mathrm{\Lambda }_s(𝜷)=\mathrm{\Lambda }[𝚽]`$ is concave up, i.e.

$$\frac{d^2\mathrm{\Lambda }_s}{d\beta _1^2}=\frac{dQ_1}{d\beta _1}>0.$$
(6)
Under this condition, the linear eigenvalue problem (4) has no unstable eigenvalues, i.e. those with a positive real part $`\lambda `$. Otherwise, the second variation (3) constrained by the set (5) has a single negative direction that corresponds to a single positive eigenvalue $`\lambda `$ in the linear eigenvalue problem (4) . The stability criterion for scalar (or one-component) NLS solitons holds when the self-adjoint operator $`𝐋_1`$ has a single negative eigenvalue, i.e. when the second variation (3), without the constraint (5) imposed, has a single negative direction. If the latter condition is not satisfied, as it happens for solitary waves with nodes, the fundamental criterion for the soliton instability can be extended only for a special case , while more generic mechanisms of oscillatory instabilities, associated with complex eigenvalues of the linear eigenvalue problem, may appear beyond the prediction of the fundamental criterion .

Here we extend the soliton stability analysis to the case of multi-component solitary waves described by a system of incoherently coupled NLS equations (1). We assume that the number of negative directions (eigenvalues) of the second variation $`\delta ^2\mathrm{\Lambda }`$ is fixed, and we denote it as $`n(\mathrm{\Lambda })`$. The unstable eigenvalues $`\lambda `$ of the linear problem (4) are connected with some negative eigenvalues of the matrix $`𝐔`$ defined by the elements

$$U_{nm}=\frac{^2\mathrm{\Lambda }_s}{\beta _n\beta _m}=\frac{Q_n}{\beta _m}=\frac{Q_m}{\beta _n}.$$
(7)
The matrix $`𝐔`$ is the Hessian matrix of the energetic surface $`\mathrm{\Lambda }_s(𝜷)`$. We denote the number of positive eigenvalues of the matrix $`𝐔`$ as $`p(U)`$, and the number of its negative eigenvalues as $`n(U)`$, so that $`p(U)+n(U)N`$, since some eigenvalues may be zeros in a degenerate (bifurcation) case. As is shown below, both $`p(U)`$ and $`n(U)`$ satisfy some additional constraints,

$$p(U)\mathrm{min}\{N,n(\mathrm{\Lambda })\},n(U)\mathrm{max}\{0,Nn(\mathrm{\Lambda })\}.$$
(8)
Within these notations, we formulate (and prove below) the following fundamental results on stability and instability of multi-component solitary waves of the coupled NLS equations (1):

(i) the linear problem (4) may have at most $`n(\mathrm{\Lambda })`$ unstable eigenvalues $`\lambda `$, all real and positive;

(ii) a multi-component soliton is linearly unstable provided $`p(U)<n(\mathrm{\Lambda })`$; then the linear problem (4) has $`n(\mathrm{\Lambda })p(U)`$ real (positive or zero-becoming-positive) eigenvalues $`\lambda `$;

(iii) a multi-component soliton is linearly stable provided $`p(U)=n(\mathrm{\Lambda })(N)`$; in the case $`n(\mathrm{\Lambda })=N`$ this criterion implies that the energetic surface $`\mathrm{\Lambda }_s(𝜷)`$ is concave up in the $`𝜷`$-space;

(iv) a single eigenvalue $`\lambda `$ crosses a marginal stability curve when the matrix $`𝐔`$ possesses a zero-becoming-negative eigenvalue; the normal form for the instability-induced dynamics of multi-component solitary waves resembles the equation of motion for an effective classical particle subjected to a $`N`$-dimensional potential field,

$$E=\frac{1}{2}\underset{n=1}{\overset{N}{}}\underset{m=1}{\overset{N}{}}M_{nm}\frac{d\nu _n}{dz}\frac{d\nu _m}{dz}+W(𝜷,𝝂),$$
(9)
where $`M_{nm}`$ are the elements of the positive-definite “mass matrix” \[see Eq. (28) below\], $`𝝂`$ is the vector describing a perturbation to the soliton parameters $`𝜷`$, and $`W(𝜷,𝝂)`$ is an effective potential energy defined as

$$W(𝜷,𝝂)=H_s(𝜷+𝝂)H_s(𝜷)+\underset{n=1}{\overset{N}{}}\left(\beta _n+\nu _n\right)\left[Q_{sn}(𝜷+𝝂)Q_{sn}(𝜷)\right].$$
(10)
These results should be compared with the results following from the stability and instability theorems earlier formulated by Grillakis, Shatah, and Strauss . The stability result (iii), i.e. the condition $`p(U)=n(\mathrm{\Lambda })`$, is identical to that of the stability theorem , but the instability results (i)-(ii) are more general and explicit. In particular, the theorem of Grillakis et al.  guarantees the soliton instability provided the difference $`n(\mathrm{\Lambda })p(U)`$ is odd. However, our results predict that the soliton instability always occurs for $`n(\mathrm{\Lambda })p(U)>0`$, being associated with exactly $`n(\mathrm{\Lambda })p(U)`$ non-negative real eigenvalues $`\lambda `$ of the linear eigenvalue problem (4). Moreover, according to our result (iv), each new unstable eigenvalue $`\lambda `$ appears via a bifurcation at the marginal stability curve where the determinant of the matrix $`𝐔`$ vanishes, i.e. it is connected with a zero-becoming-negative eigenvalue of the Hessian matrix $`𝐔`$. If $`n(\mathrm{\Lambda })>N`$, unstable eigenvalues $`\lambda `$ originated from the negative eigenvalues of the Hessian matrix $`𝐔`$ co-exist with $`n(\mathrm{\Lambda })N`$ unstable eigenvalues of the linear problem (4), i.e. a solitary wave is unconditionally unstable when $`n(\mathrm{\Lambda })>N`$.

## III A Proof of the Basic Results

Here we develop the analysis of the problem (3)-(5), in order to prove the results (i)-(iv) formulated above. The Sturm-Liouville operators $`(L_0)_{nn}`$ are all non-negative since the fundamental (node-less) solution $`𝚽(𝐱)`$ for a solitary wave reaches a bottom of the spectrum at zero: $`(L_0)_{nn}\mathrm{\Phi }_n=0`$. As a result,

$$\text{(min) }\delta ^2\mathrm{\Lambda }=_{\mathrm{}}^{\mathrm{}}𝑑𝐱𝐮|𝐋_1𝐮=\underset{\mu }{}\mu _{\mathrm{}}^{\mathrm{}}𝑑𝐱𝐮_k|𝐮_k.$$
(11)
Here ($`\mu `$,$`𝐮_k`$) are eigenvalues and eigenfunctions of the auxiliary linear problem,

$$𝐋_1𝐮_k=\mu 𝐮_k\underset{m=1}{\overset{N}{}}\nu _m\mathrm{\Phi }_m(𝐱)𝐞_m.$$
(12)
The linear problem (12) is constrained by the set (5) and the parameters $`\nu _1`$, $`\nu _2`$,…,$`\nu _N`$ have the meaning of the Lagrange multipliers.

Let us suppose that the Sturm-Liouville matrix operator $`𝐋_1`$ has $`n(\mathrm{\Lambda })`$ negative eigenvalues $`\mu =\{\mu _{n(\mathrm{\Lambda })},\mu _{n(\mathrm{\Lambda })+1},\mathrm{},\mu _1\}`$ corresponding to the eigenfunctions $`𝐮=\{𝝍_{n(\mathrm{\Lambda })}(𝐱),𝝍_{n(\mathrm{\Lambda })+1}(𝐱),\mathrm{},𝝍_1(𝐱)\}`$; a single zero eigenvalue with a one-node eigenfunction $`𝐮=d𝚽/dx`$ and the rest of the spectrum is positive and contains $`N`$ branches of the continuous spectrum, for $`\mu >\{\beta _1,\beta _2,\mathrm{},\beta _N\}`$, and some isolated positive eigenvalues, for $`\mu =\{\mu _1,\mu _2,\mathrm{},\mu _p\}`$.

The mathematical problem can then be reformulated in the following way. The linear operator $`𝐋_1`$ has $`n(\mathrm{\Lambda })`$ negative eigenvalues that generate negative directions of the second variation $`\delta ^2\mathrm{\Lambda }`$. However, the corresponding eigenfunctions do not satisfy generally the constraints (5). By introducing the Lagrangian multipliers in (11) and (12), we satisfy a constrained minimization problem (3) and (5) but, due to this procedure, the number of negative eigenvalues may be reduced. We will show how to connect the total number of negative eigenvalues of the constrained problem (5), (11), and (12) with the negative eigenvalues of the Hessian matrix (7). But, as a prerequisite, we prove two additional results for the spectrum of the problem (4):

(i) the spectrum of $`\lambda ^2`$ is real, i.e. oscillatory instabilities are prohibited;

(ii) each negative direction ($`\mu ,𝐮_k`$) of the problem (12) generates an unstable (positive) eigenvalue $`\lambda `$ of the problem (4).

To prove the statement (i), we notice that the matrix operator $`𝐋_0`$ can be factorized as $`𝐋_0=_{d=1}^D𝐌_d^+𝐌_d`$, where $`𝐌_d`$ has a diagonal form with the following matrix elements

$$\left(M_d\right)_{nn}=\sqrt{d_n}\left[_{x_d}+\frac{1}{\mathrm{\Phi }_n(𝐱)}_{x_d}\mathrm{\Phi }_n(𝐱)\right],$$

provided the soliton solutions $`\mathrm{\Phi }_n(𝐱)`$ have no nodes in a finite domain. Using this factorization, the linear problem (4) can be rewritten for the function $`𝐮=_{d=1}^D𝐌_d^+𝐯_d`$ as follows,

$$\underset{d^{}=1}{\overset{D}{}}𝐌_d𝐋_1𝐌_d^{}^+𝐯_d^{}=\lambda ^2𝐯_d.$$

Since the matrix operator with the elements $`𝐌_d𝐋_1𝐌_d^{}^+`$ is Hermitian, its eigenvalues ($`\lambda ^2`$) are all real.

To prove the statement (ii), we suppose that we have constructed a negative direction ($`\mu ,𝐮_k`$) of the problem (11)-(12) subject to the constraints (5). Then, the linear problem (4) has an unstable eigenvalue $`\lambda `$ defined as,

$$\lambda ^2=\frac{𝐮_k|𝐋_0𝐋_1𝐮_k}{𝐮_k|𝐮_k}=\mu \frac{𝐮_k|𝐋_0𝐮_k}{𝐮_k|𝐮_k}.$$
(13)
Since the linear operator $`𝐋_0`$ is positive definite for any $`𝐮_k(\mathrm{𝟎},𝚽)`$, we have $`\lambda ^2>0`$ for any $`\mu <0`$.

Our next goal is to construct solutions to the auxillary problem (12). Since the matrix operator $`𝐋_1`$ is Hermitian, it has a complete spectrum in a Hilbert space that is suitable for expanding the vector function $`𝐮_k(𝐱)`$. We present such a spectral decomposition in the form,

$$𝐮_k(𝐱)=\underset{m=1}{\overset{N}{}}\nu _m\left(\underset{\mu _r<0}{}\frac{𝝍_r|\mathrm{\Phi }_m𝐞_m}{\mu \mu _r}𝝍_r(𝐱)+\underset{\mu _r>0}{}\frac{𝝍_r|\mathrm{\Phi }_m𝐞_m}{\mu \mu _r}𝝍_r(𝐱)\right),$$
(14)
where the sum $`_{\mu _r<0}`$ contains $`n(\mathrm{\Lambda })`$ terms from the negative spectrum, while the sum $`_{\mu _r>0}`$ includes schematically both the discrete and continuous positive spectra of $`𝐋_1`$. The contribution from the neutral eigenfunction $`𝐮=d𝚽/dx`$ vanishes due to the symmetry properties. General solution (14) has to be constrained by the conditions (5). This system reduces to the linear algebra for the Lagrange multipliers, $`𝐀(\mu )\nu =\mathrm{𝟎}`$, where the matrix $`𝐀(\mu )`$ has a symmetric form with the elements:

$$A_{nm}(\mu )=\underset{\mu _r<0}{}\frac{\mathrm{\Phi }_n𝐞_n|𝝍_r𝝍_r|\mathrm{\Phi }_m𝐞_m}{\mu \mu _r}+\underset{\mu _r>0}{}\frac{\mathrm{\Phi }_n𝐞_n|𝝍_r𝝍_r|\mathrm{\Phi }_m𝐞_m}{\mu \mu _r}.$$
(15)
The linear system $`𝐀(\mu )𝝂=\gamma 𝝂`$ has generally $`N`$ real eigenvalues $`\gamma _1(\mu )`$, $`\gamma _2(\mu )`$, …, $`\gamma _N(\mu )`$. These eigenvalues are continuous functions of $`\mu `$ for $`\mu 0`$, except for $`n(\mathrm{\Lambda })`$ resonant planes at $`\mu =\{\mu _{n(\mathrm{\Lambda })},\mu _{n(\mathrm{\Lambda })+1},\mathrm{},\mu _1\}`$. At these planes, the matrix $`𝐀(\mu )`$ has poles and the eigenvalues $`\gamma (\mu )`$ may have singularities. Below, we prove the following three properties of the eigenvalues $`\gamma (\mu )`$:

(i) all eigenvalues $`\gamma (\mu )`$ are negative for $`\mu <\mu _{n(\mathrm{\Lambda })}(<0)`$;

(ii) each eigenvalue $`\gamma (\mu )`$ is a decreasing function of $`\mu `$ for $`\mu 0`$, except for $`n(\mathrm{\Lambda })`$ resonant planes at $`\mu =\{\mu _{n(\mathrm{\Lambda })},\mu _{n(\mathrm{\Lambda })+1},\mathrm{},\mu _1\}`$;

(iii) at least ($`N1`$) eigenvalues $`\gamma (\mu )`$ are continuous at any of the resonant planes $`\mu =\mu _r<0`$, while the minimal eigenvalue, say $`\gamma _1(\mu )`$, may have an infinite discontinuity jumping from negative infinity, at $`\mu \mu _r0`$, to positive infinity, at $`\mu \mu _r+0`$.

To show the property (i), we consider the asymptotic limit of $`𝐀(\mu )`$ as $`\mu \mathrm{}`$. In this limit, the eigenvalues $`\gamma (\mu )`$ can be expressed from the algebra of quadratic forms,

$$\gamma (\mu )=\frac{1}{\mu 𝝂|𝝂}\left(\underset{\mu _r<0}{}b_r+\underset{\mu _r>0}{}b_r\right),$$
(16)
where

$$b_r=\left|\underset{n=1}{\overset{N}{}}\nu _n𝝍_r|\mathrm{\Phi }_n𝐞_n\right|^20.$$
(17)
Since all $`b_r`$ may not vanish simultaneously for $`𝝂\mathrm{𝟎}`$, the eigenvalues $`\gamma (\mu )`$ are negative definite in Eq. (16) so that $`\gamma (\mu )0`$ as $`\mu \mathrm{}`$.

To show the property (ii), we take the derivative of the system $`𝐀𝝂=\gamma (\mu )𝝂`$ and use the algebra of quadratic forms. The derivative of $`\gamma (\mu )`$ is then defined for $`\mu 0`$ excluding the resonant planes at $`\mu =\{\mu _{n(\mathrm{\Lambda })},\mu _{n(\mathrm{\Lambda })+1},\mathrm{},\mu _1\}`$ as

$$\frac{d\gamma (\mu )}{d\mu }=\frac{1}{𝝂|𝝂}𝝂|\frac{d𝐀(\mu )}{d\mu }𝝂=\frac{1}{𝝂|𝝂}\left(\underset{\mu _r<0}{}\frac{b_r}{(\mu \mu _r)^2}+\underset{\mu _r>0}{}\frac{b_r}{(\mu \mu _r)^2}\right),$$
(18)
where $`b_r`$ are defined by the same relation (17). Since the derivative of $`\gamma (\mu )`$ is hegative definite in Eq. (18), all eigenvalues $`\gamma (\mu )`$ are decreasing functions of $`\mu `$ whenever $`d\gamma (\mu )/d\mu `$ exists.

To show the property (iii), we consider the behavior of the eigenvalues $`\gamma (\mu )`$ at the resonant plane $`\mu =\mu _r<0`$. In this limit, the matrix elements $`A_{nm}(\mu )`$ have the following asymptotic form,

$$A_{nm}(\mu )\frac{\mathrm{\Phi }_n𝐞_n|𝝍_r𝝍_r|\mathrm{\Phi }_m𝐞_m}{(\mu \mu _r)}.$$

Therefore, the matrix $`𝐀(\mu )`$ has ($`N1`$) zero eigenvalues $`\gamma (\mu )`$ and a single non-zero eigenvalue $`\gamma _1(\mu )`$ with the asymptotic value,

$$\gamma _1(\mu )\frac{1}{(\mu \mu _r)}\underset{n=1}{\overset{N}{}}\left|\mathrm{\Phi }_n𝐞_n|𝝍_r\right|^2.$$
(19)
If the sum in Eq. (19) does not vanish, the eigenvalue $`\gamma _1(\mu )`$ has an infinite discontinuity described in (iii) and, according to the property (ii), it is the minimal eigenvalue. Other ($`N1`$) eigenvalues are in fact non-zero in the limit $`\mu \mu _r`$. Since the matrix $`𝐀(\mu )`$ is a meromorphic function of $`\mu `$ as $`\mu 0`$, the eigenvalue $`(\mu \mu _r)\gamma (\mu )`$ is of order of $`O(\mu \mu _r)`$ for ($`N1`$) non-singular eigenvalues. Therefore, the values of $`\gamma (\mu )`$ are generally non-zero in the limit $`\mu \mu _r`$.

Thus, we have a clear picture how the eigenvalues $`\gamma (\mu )`$ behave as functions of $`\mu `$ \[see Figs 1(a,b)\]. Starting with small negative values as $`\mu \mathrm{}`$, all eigenvalues decrease as $`\mu `$ grows towards the $`n(\mathrm{\Lambda })`$ resonant planes. At each of those planes, ($`N1`$) eigenvalues remain continuously decreasing, while one (minimal) eigenvalue jumps to a positive domain unless the condition

$$\underset{n=1}{\overset{N}{}}\left|\mathrm{\Phi }_n𝐞_n|𝝍_r\right|^2=0$$
(20)
is satisfied (this condition will be discussed below). Assuming that the condition (20) is not met, we come to the conclusion that a root of $`\gamma (\mu )`$ may occur only after a jump of $`\gamma (\mu )`$ at a resonant plane $`\mu =\mu _r`$ to a large positive value, and further decrease of $`\gamma (\mu )`$ as $`\mu (>\mu _r)`$ grows. The root of $`\gamma (\mu )`$, if exists for $`\mu 0`$, produces a legitimate solution $`𝐮_k(𝐱)`$ of the problem (12) under the constraints (5). This solution ($`\mu ,𝐮_k`$) would then be associated with an unstable eigenvalue $`\lambda `$, according to the connection formula (13). Thus, our main task is to control the behavior of positive $`\gamma (\mu )`$ in between the plane $`\mu =0`$ and the resonant planes $`\mu =\{\mu _{n(\mathrm{\Lambda })},\mu _{n(\mathrm{\Lambda })+1},\mathrm{},\mu _1\}`$.

At the plane $`\mu =0`$, the problem (12) has a simple solution for $`𝐮_k(𝐱)`$,

$$𝐮_{\mu =0}(𝐱)=\underset{n=1}{\overset{N}{}}\nu _n\frac{𝚽(𝐱)}{\beta _n}.$$
(21)
Substituting Eq. (21) into the constraints (5), we find that $`𝐀(0)=𝐔`$, where the matrix $`𝐔`$ is the Hessian of the energetic surface $`\mathrm{\Lambda }_s(\beta )`$ with the elements $`U_{nm}`$ defined by Eq. (7). We can now use this construction and prove the main results (i)–(iii). In the analysis below we assume that the condition (20) is never met and the root of $`\gamma (\mu )`$ at $`\mu 0`$ is associated with the unstable eigenvalue $`\lambda `$ of the stability problem (4).

The roots of $`\gamma (\mu )`$ may appear only to the right of any of the $`n(\mathrm{\Lambda })`$ resonant plane. There are totally $`n(\mathrm{\Lambda })`$ jumps of $`\gamma (\mu )`$ to positive values at $`\mu 0`$ and, therefore, no more than $`n(\mathrm{\Lambda })`$ roots of $`\gamma (\mu )`$ may exist for $`\mu 0`$.

If $`n(\mathrm{\Lambda })=N`$, the positive eigenvalues $`\gamma (\mu )`$ remain continuous after passing the corresponding resonant plane at $`\mu =\mu _r(<0)`$. Therefore, the sign of these eigenvalues is controlled by the eigenvalues of the Hessian matrix $`𝐔`$ at $`\mu =0`$. If $`p(U)=N=n(\mathrm{\Lambda })`$, all positive eigenvalues $`\gamma (\mu )`$ remain positive for $`\mu _r<\mu 0`$ and no roots of $`\gamma (\mu )`$ exist for $`\mu 0`$ \[see Fig. 1(a)\]. If $`p(U)<N=n(\mathrm{\Lambda })`$, there exist $`Np(U)`$ negative or zero-becoming-negative eigenvalues of $`𝐔`$ that correspond to $`Np(U)`$ roots of $`\gamma (\mu )`$ for $`\mu 0`$.

If $`n(\mathrm{\Lambda })<N`$, then $`Nn(\mathrm{\Lambda })`$ eigenvalues $`\gamma (\mu )`$ do not have jumps at the corresponding resonant planes $`\mu =\mu _r`$. They continue to be negative and match at $`\mu =0`$ with the $`Nn(\mathrm{\Lambda })`$ negative eigenvalues of $`𝐔`$. From this, we come to the conclusion that $`p(U)`$ and $`n(U)`$ satisfy the constraints (8), i.e. $`n(U)Nn(\mathrm{\Lambda })`$ or, equivalently, $`p(U)n(\mathrm{\Lambda })`$. Furthermore, other $`n(\mathrm{\Lambda })`$ $`(<N)`$ eigenvalues $`\gamma (\mu )`$ may have roots for $`\mu 0`$ which are completely controlled by the remaining $`n(\mathrm{\Lambda })`$ eigenvalues of $`𝐔`$ according to the same criterion as in the case $`p(\mathrm{\Lambda })=N`$. For instance, if $`p(U)<n(\mathrm{\Lambda })`$, then $`n(\mathrm{\Lambda })p(U)`$ negative or zero-becoming-negative eigenvalues of the matrix $`𝐔`$ correspond to $`n(\mathrm{\Lambda })p(U)`$ roots of $`\gamma (\mu )`$ at $`\mu 0`$.

If $`n(\mathrm{\Lambda })>N`$, then $`n(\mathrm{\Lambda })N`$ eigenvalues $`\gamma (\mu )`$ jump twice in the domain $`\mu 0`$ leading to at least $`n(\mathrm{\Lambda })N`$ unconditional roots for $`\mu 0`$ \[see Fig.1(b)\]. After the jumps, the $`N`$ eigenvalues $`\gamma (\mu )`$ match the $`N`$ eigenvalues of the matrix $`𝐔`$ and may have additional roots of $`\gamma (\mu )`$ if $`p(U)<N`$. Total number of roots of $`\gamma (\mu )`$ at $`\mu 0`$ is then defined as $`(n(\mathrm{\Lambda })N)+(Np(U))=n(\mathrm{\Lambda })p(U)`$.

The analysis above is valid for the non-degenerate case when the condition (20) is never satisfied. However, the stability and instability results (i)-(iii) are not affected even if the condition (20) is satisfied for a particular resonant plane $`\mu =\mu _r(<0)`$. In this case, the eigenfunction $`𝐮_k(𝐱)`$ of the operator $`𝐋_1`$ satisfies all the constraints (5) identically and, therefore, the eigenvalue $`\mu =\mu _r`$ is associated with an unstable eigenvalue $`\lambda `$, according to Eq. (13). Although the eigenvalue $`\gamma _1(\mu )`$ has no jump at $`\mu =\mu _r`$ \[see Eq. (19)\] and is continuous, it is still controlled by the negative eigenvalues of $`𝐔`$ at $`\mu =0`$. Indeed, in this case, the minimal eigenvalue $`\gamma _1(\mu )`$ at $`\mu <\mu _r`$ remains negative for $`\mu >\mu _r`$ and matches with a negative eigenvalue of the matrix $`𝐔`$ (if no other jumps occur in the domain $`\mu 0`$). This additional negative eigenvalue $`\mu `$ still predicts the instability, according to the result (iii).

Finally, we prove the result (iv) for the instability bifurcation of multi-component solitary waves. Provided the number $`n(\mathrm{\Lambda })`$ is fixed, the instability bifurcation may only occur when $`𝐀(0)=𝐔`$ has a zero eigenvalue for a certain eigenvector $`𝝂=𝝂^{(k)}`$. Let us define $`𝐔=𝐔_{\mathrm{thr}}`$ at the marginal stability curve so that the determinant of $`𝐔_{thr}`$ vanishes. The instability bifurcations of multi-component solitons were considered in Refs.  but the results do not agree with each other. Here, we recover the results of Ref.  and derive the normal form (9) by an elegant reduction of general algebraic expressions.

Assuming $`\mu =0`$ for $`𝝂=𝝂^{(k)}`$ so that $`𝐔_{thr}𝝂^{(k)}=\mathrm{𝟎}`$, we find the asymptotic solution of Eq. (4) in the form (21) and

$$𝐰_{\mu =0}(𝐱)=\lambda \underset{n=1}{\overset{N}{}}\nu _n^{(k)}𝐋_0^1\frac{𝚽(𝐱)}{\beta _n}.$$
(22)
In this limit, the second variation $`\delta ^2\mathrm{\Lambda }`$ of the Lyapunov functional can be found from Eqs. (3), (5), (21), and (22) as follows

$$\delta ^2\mathrm{\Lambda }=D_1\lambda ^2,$$
(23)
where

$$D_1=_{\mathrm{}}^{\mathrm{}}𝑑𝐱\underset{m=1}{\overset{N}{}}\frac{1}{\mathrm{\Phi }_m^2(𝐱)}\left[\underset{n=1}{\overset{N}{}}\nu _n^{(k)}_0^𝐱𝑑𝐱^{}\mathrm{\Phi }_m(𝐱^{})\frac{\mathrm{\Phi }_m(𝐱^{})}{\beta _n}\right]^2.$$
(24)
The integral converges under the condition that $`𝝂^{(k)}`$ is a solution of the equation $`𝐔_{\mathrm{thr}}𝝂^{(k)}=\mathrm{𝟎}`$. On the other hand, the perturbation (21) shifts the soliton parameter $`𝜷`$ according to the expression, $`𝚽(𝐱,𝜷)+𝐮_{\mu =0}(𝐱)𝚽(𝐱,𝜷+𝝂^{(k)})`$. As a result, the second variation can be closed as

$$\delta ^2\mathrm{\Lambda }=2\left[\mathrm{\Lambda }\mathrm{\Lambda }_s(𝜷+𝝂^{(k)})\right]D_0,$$
(25)
where

$$D_0=𝝂^{(k)}|𝐔𝝂^{(k)}.$$
(26)
The parameter $`\mathrm{\Lambda }`$ in Eq. (25) is chosen from the condition that the first variation of $`\mathrm{\Lambda }_s(𝜷+𝝂^{(k)})`$ vanishes for arbitrary $`𝝂^{(k)}`$. This gives the connection formula: $`\mathrm{\Lambda }\mathrm{\Lambda }_{st}=H_s(𝜷)+_{n=1}^N(\beta _n+\nu _n^{(k)})Q_{sn}(𝜷)`$. Equating Eq. (23) and Eq. (25), we recover the result of the bifurcation theory,

$$\lambda ^2=\frac{D_0}{D_1}.$$
(27)
Since $`D_1>0`$ \[see Eq. (24)\], the positive values of $`\lambda ^2`$ occur when the determinant of the matrix $`𝐔`$ is small and negative (i.e. the matrix $`𝐔`$ has a zero-becoming-negative eigenvalue when the soliton parameter $`𝜷`$ crosses the marginal stability curve). The explicit formulas of the soliton bifurcation theory provide an alternative and more compact form for the determinants $`D_0`$ and $`D_1`$ compared to those obtained in Ref. .

The normal form (9) follows from Eqs. (23) and (25) when $`\mathrm{\Lambda }=\mathrm{\Lambda }_{st}+E`$, and the perturbation vector $`𝝂^{(k)}`$ is replaced by a slowly varying vector $`𝝂=𝝂(z)`$ (see  for details). Then, the surface $`\mathrm{\Lambda }_s(𝜷+𝝂)`$ is extended beyond the second variation limit, and the linear approximation is converted into the slope: $`\lambda 𝝂(z)=d𝝂(z)/dz`$. The mass constants $`M_{nm}`$ follow from Eq. (24) in the explicit form:

$$M_{nm}=_{\mathrm{}}^{\mathrm{}}𝑑𝐱\underset{k=1}{\overset{N}{}}\frac{1}{\mathrm{\Phi }_k^2(𝐱)}\left(_0^𝐱𝑑𝐱^{}\mathrm{\Phi }_k(𝐱^{})\frac{\mathrm{\Phi }_k(𝐱^{})}{\beta _n}\right)\left(_0^𝐱𝑑𝐱^{}\mathrm{\Phi }_k(𝐱^{})\frac{\mathrm{\Phi }_k(𝐱^{})}{\beta _m}\right).$$
(28)
The normal form (9) resembles the conserved sum of the kinetic energy and potential energy $`W(𝜷,𝝂)`$ of a particle moving in a $`N`$-dimensional space. We notice that the kinetic energy with the “mass” matrix (28) is positive definite and the unperturbed multicomponent solitary wave (i.e. that with $`𝝂=\mathrm{𝟎}`$) is a stationary point of $`W(𝜷,𝝂)`$ for any $`𝝂`$. Thus, the stability of multi-component solitons resembles the stability of a particle located at an equilibrium point of the $`N`$-dimensional field . Under the condition that $`n(\mathrm{\Lambda })=N`$, the particle is stable if in the $`𝜷`$-space the potential energy surface $`W(𝜷)`$ is concave up, and it is unstable if the potential energy surface is saddle-type or concave down. If $`n(\mathrm{\Lambda })<N`$, the potential energy surface $`W(𝜷)`$ always has some $`Nn(\mathrm{\Lambda })`$ negative directions that do not affect the stability properties of the particle. However, the remaining $`n(\mathrm{\Lambda })`$ $`(<N)`$ directions of the potential energy surface define the stability of the particle with the same criterion as above. Finally, for the case $`n(\mathrm{\Lambda })>N`$, the soliton stability properties defined by the type of the potential energy surface $`W(𝜷)`$ are not conclusive since the corresponding unstable eigenvalues coexist with additional $`n(\mathrm{\Lambda })N`$ unconditionally unstable eigenvalues.

## IV Example: two coupled NLS equations

In order to demonstrate how our general theory can be applied to a particular physical problem and also to compare the stability and instability results (ii)-(iii) with some earlier known examples, we consider here an important case of two incoherently coupled NLS equations in (1+1) dimension (see, e.g., ):

$`i{\displaystyle \frac{\psi _1}{z}}+{\displaystyle \frac{^2\psi _1}{x^2}}+\left(|\psi _1|^2+\gamma |\psi _2|^2\right)\psi _1`$ $`=`$ $`0,`$ (29)
$`i{\displaystyle \frac{\psi _2}{z}}+{\displaystyle \frac{^2\psi _2}{x^2}}+\left(|\psi _2|^2+\gamma |\psi _1|^2\right)\psi _2`$ $`=`$ $`0,`$ (30)
where $`\gamma `$ is a coupling parameter. System (30) is a two-component reduction of the general $`N`$component system (1) for $`d_1=d_2=1`$, $`\gamma _{11}=\gamma _{22}=1`$, and $`\gamma _{12}=\gamma _{21}=\gamma `$. An explicit soliton solution can be easily found for $`\beta _1=\beta _2=\beta `$ and $`\gamma >1`$ in the form,

$$\mathrm{\Phi }_1(x)=\mathrm{\Phi }_2(x)=\sqrt{\frac{2\beta }{1+\gamma }}\mathrm{sech}\left(\sqrt{\beta }x\right).$$
(31)
This solution describes a two-component solitary wave with the components of equal amplitudes. It corresponds to a straight line $`\beta _1=\beta _2`$ in the parameter plane ($`\beta _1,\beta _2`$) of a general two-parameter family of solitary waves of the model (30). When $`1<\gamma 0`$, such two-parameter solitons may exist everywhere in the plane ($`\beta _1,\beta _2`$), while for $`\gamma >0`$, the soliton existence domain is restricted by two bifurcation curves, $`\beta _2=\omega _\pm (\gamma )\beta _1`$, where

$$\omega _\pm (\gamma )=\left(\frac{\sqrt{1+8\gamma }1}{2}\right)^{\pm 2}.$$
(32)
Approximate analytical expressions can also be obtained in the vicinity of the bifurcation curves (32), when one of the components of a composite solitary wave becomes small, while the other one is described by the scalar NLS equation. Such a case, when one of the component creates an effective waveguide that guides the other component, is known to describe the so-called shepherding effect where the large-amplitude component plays a role of a shepherding pulse . The composite soliton, that describes a shepherding pulse $`\psi _1`$ guiding a small component pulse $`\psi _2`$, can be found in the form (see also Ref. ),

$$\mathrm{\Phi }_1=R_0(x)+ϵ^2R_2(x)+\mathrm{O}(ϵ^4),\mathrm{\Phi }_2=ϵS_1(x)+\mathrm{O}(ϵ^3).$$
(33)
It exists in the vicinity of the bifurcation curve,

$$\beta _2=\omega _+(\gamma )\beta _1+ϵ^2\omega _{2+}(\gamma )\beta _1+\mathrm{O}(ϵ^4),$$
(34)
and the main terms of the asymptotic series (33), (34) are defined as

$$R_0=\sqrt{2\beta _1}\mathrm{sech}\left(\sqrt{\beta _1}x\right),S_1=\sqrt{\beta _1}\mathrm{sech}^{\sqrt{\omega _+}}\left(\sqrt{\beta _1}x\right),$$

and

$$\omega _{2+}=\frac{_{\mathrm{}}^{\mathrm{}}𝑑x\left(S_1^4+2\gamma R_0R_2S_1^2\right)}{_{\mathrm{}}^{\mathrm{}}𝑑xS_1^2},$$

The second-order correction $`R_2(x)`$ is a solution of the differential equation,

$$\left[_x^2+\beta _16\beta _1\mathrm{sech}^2\left(\sqrt{\beta _1}x\right)\right]R_2=\gamma R_0S_1^2.$$

From the domain of existence of the two-component soliton, it follows that $`\omega _{2+}(\gamma )>0`$, for $`0<\gamma <1`$, and $`\omega _{2+}(\gamma )<0`$, for $`\gamma >1`$. At $`\gamma =1`$ (the so-called integrable Manakov case), a family of two-parameter composite solitons becomes degenerated: it exists on the line $`\beta _1=\beta _2`$ but, generally, it is different from the one-parameter solution (31). The coupled solitons are known to be stable for the integrable case $`\gamma =1`$. Here we apply the stability theory developed above and prove that the (1+1)-dimensional two-parameter solitons, including the solitons of equal amplitudes (31), are stable for $`\gamma 0`$, and unstable for $`\gamma <0`$.

First, we evaluate the indices $`p(U)`$ and $`n(\mathrm{\Lambda })`$ for the explicit solution (31). As follows from Eqs. (30) and (31), the Hessian matrix $`𝐔`$ with the elements (7) can be found in the form,

$$\frac{Q_1}{\beta _1}=\frac{Q_2}{\beta _2}=\frac{1}{\sqrt{\beta }(1+\gamma )}\text{and}\frac{Q_1}{\beta _2}=\frac{Q_2}{\beta _1}=\frac{\gamma }{\sqrt{\beta }(1+\gamma )}.$$

As follows from these results, the Hessian matrix has $`p(U)=2`$ positive eigenvalues, for $`1<\gamma <1`$, and $`p(U)=1`$ positive eigenvalue, for $`\gamma >1`$. On the other hand, the linear matrix operator $`𝐋_1`$ given below Eq. (3) can be diagonalized for linear combinations of the eigenfunctions, $`v_1=u_1+u_2`$ and $`v_2=u_1u_2`$, such that

$`\left[_x^2+\beta 6\beta \mathrm{sech}^2\left(\sqrt{\beta }x\right)\right]v_1`$ $`=`$ $`\mu v_1,`$ (35)
$`\left[_x^2+\beta 2\beta {\displaystyle \frac{(3\gamma )}{(1+\gamma )}}\mathrm{sech}^2\left(\sqrt{\beta }x\right)\right]v_2`$ $`=`$ $`\mu v_2.`$ (36)
Both the operators in Eqs. (36) are the linear Schrödinger operators with solvable sech-type potentials, and the corresponding eigenvalue spectra are well studied. The first operator always has a single negative eigenvalue for $`\mu =3\beta `$, whereas the second operator has no negative eigenvalues, for $`\gamma >1`$, has a single negative eigenvalue, for $`0<\gamma <1`$, and it has two negative eigenvalues, for $`1<\gamma <0`$. Thus, in total there exist $`n(\mathrm{\Lambda })=3`$ negative eigenvalues, for $`1<\gamma <0`$, $`n(\mathrm{\Lambda })=2`$ negative eigenvalues, for $`0<\gamma <1`$, and $`n(\mathrm{\Lambda })=1`$ negative eigenvalue, for $`\gamma >1`$.

Applying the stability and instability results (ii)-(iii) obtained and discussed in Secs. II and III, we come to the conclusion that the soliton solution (31) with equal amplitudes is linearly stable for $`\gamma >0`$, since in this domain $`p(U)=n(\mathrm{\Lambda })=\{1,2\}`$, and it is linearly unstable for $`1<\gamma <0`$, since in this domain $`p(U)=2<n(\mathrm{\Lambda })=3`$.

The soliton stability in the model (30) for $`\gamma >0`$ was also studied by Bergé  who considered the case of degenerate one-parametric solitary waves (31). Here we have extended those results to a general case: the same stability and instability results are valid for the two-parameter family of solitons provided that the indices $`p(U)`$ and $`n(\mathrm{\Lambda })`$ remain unchanged for the values ($`\beta _1,\beta _2`$) that do not belong to the line $`\beta _1=\beta _2`$. Indeed, applying a perturbation theory for small $`\gamma `$ (see  for details), one can show that $`n(\mathrm{\Lambda })=3`$, for $`1\gamma <0`$, and $`n(\mathrm{\Lambda })=2`$, for $`0<\gamma 1`$.

Finally, we consider the other limiting case that describes the shepherding effect \[see Eqs. (33) and (34)\]. In this limit, the elements (7) of the Hessian matrix $`𝐔`$ can also be calculated in explicit analytical form,

$$\frac{Q_1}{\beta _1}=\frac{1}{\sqrt{\beta _1}}+\frac{r_2^2}{\omega _{2+}s_1}+\mathrm{O}(ϵ^2),\frac{Q_1}{\beta _2}=\frac{Q_2}{\beta _1}=\frac{r_2}{\omega _{2+}}+\mathrm{O}(ϵ^2),\frac{Q_2}{\beta _2}=\frac{s_1}{\omega _{2+}},$$

where

$$s_1=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}S_1^2𝑑x\text{ and }r_2=_{\mathrm{}}^{\mathrm{}}R_0R_2𝑑x.$$

Since $`s_1>0`$ for any $`\gamma `$, while $`\omega _{2+}(\gamma )>0`$ for $`0<\gamma <1`$ and $`\omega _{2+}(\gamma )<0`$ for $`\gamma >1`$, the Hessian matrix $`𝐔`$ calculated for the shepherding soliton (33) has $`p(U)=2`$ positive eigenvalues, for $`0<\gamma <1`$, and $`p(U)=1`$ positive eigenvalue, for $`\gamma >1`$.

On the other hand, the linear matrix operator $`𝐋_1`$ cannot be diagonalized for the shepherding soliton (33) unless $`ϵ=0`$. In the latter (decoupled) case, it has a single negative eigenvalue at $`\mu =3\beta _1`$ and a double degenerate zero eigenvalue. When $`ϵ0`$, the zero eigenvalue shifts to become $`\mu =2\omega _{2+}(\gamma )\beta _1ϵ^2+\mathrm{O}(ϵ^4)`$. Therefore, the matrix operator $`𝐋_1`$ for the shepherding soliton (33) has $`n(\mathrm{\Lambda })=2`$ negative eigenvalues, for $`0<\gamma <1`$, and $`n(\mathrm{\Lambda })=1`$ negative eigenvalue, for $`\gamma >1`$. Thus, we come to the conclusion that the shepherding soliton is stable for $`\gamma >0`$ since $`p(U)=n(\mathrm{\Lambda })=\{1,2\}`$.

## V Conclusion

We have developed a rigorous stability analysis of multi-component solitary waves considering a system of incoherently coupled NLS equations (1) as a particular but important physical example. The method and, correspondingly, both stability and instability results can be extended to other types of solitary waves, such as multi-component spatial solitons (e.g., incoherent solitons) in non-Kerr (e.g., saturable) media, parametric solitary waves in quadratic (or $`\chi ^{(2)}`$) optical media, etc. In all such cases, our stability and instability results (i)–(iv) can be readily generalized with a rigorous proof of some of the previously known results of the asymptotic multi-scale expansion theory. However, additional analysis is required in each of those cases in order to clarify the conditions when these results completely define the stability properties of multi-component solitary waves. In the cases beyond these conditions, oscillatory instabilities may occur, and appropriate studies should rely solely on the numerical analysis of the corresponding eigenvalue problems and their linear spectra.

## Acknowledgments

The authors appreciate collaboration with A.A. Sukhorukov and E.A. Ostrovskaya. Dmitry Pelinovsky thanks C.K.R.T. Jones, K. Promislow, M. Weinstein, J. Yang, and A. Yew for stimulating discussions at different stages of the work. The work is supported by the Performance and Planning Fund of the Institute of Advanced Studies and by the Australian Photonics Cooperative Research Centre.

Figure Captions

Fig. 1: Eigenvalues $`\gamma `$ versus $`\mu `$ in the problem $`𝐀(\mu )𝝂=\gamma (\mu )𝝂`$ for $`N=3`$: (a) a stable problem with no roots of $`\gamma (\mu )`$ for $`\mu 0`$, when $`p(U)=n(\mathrm{\Lambda })=3`$; (b) an unstable problem with a single root of $`\gamma (\mu )`$ for $`\mu 0`$, when $`p(U)=3<n(\mathrm{\Lambda })=4`$.


text:
# University of Wisconsin - Madison MADPH-00-1174May 2000 High Energy Neutrinos from Gamma Ray Bursts: Event Rates in Neutrino Telescopes

## 1 Introduction

The origin of GRB is one of the most fascinating outstanding problems in astronomy. Their observed energy injection in the Universe is sufficiently large to possibly resolve another long-standing puzzle: the origin of the highest energy cosmic rays . Mounting evidence suggests that GRB emission is produced by a relativistically expanding fireball, energized by a process involving neutron stars or black holes . In the early stages the fireball, its radiation trapped by the very large optical depth, cannot emit photons efficiently. The fireball’s kinetic energy is therefore dissipated until it becomes optically thin — a scenario that can explain the observed energy and time scales of GRBs, provided the bulk Lorentz factor of the expanding flow, $`\mathrm{\Gamma }`$, is $`10^210^3`$.

Protons accelerated in shocks in the expanding fireball interact with photons to produce charged pions, the parents of high-energy neutrinos . Assuming that particles accelerated in the GRB sources produce the observed cosmic rays above the “ankle” of the energy spectrum near $`3\times 10^{18}`$ eV, one derives that the average single burst produces only $`10^2`$ neutrino events in a high energy neutrino telescope with $`1\mathrm{km}^2`$ effective area. Although the expected rate is therefore low, the neutrino signal should be relatively easy to observe provided the detector is large enough: GRB neutrinos will have a hard spectrum extending well beyond the background from atmospheric neutrinos and, even more important, the high-energy GRB neutrino events should coincide with observed GRB photon events within a narrow time window.

In this paper we calculate the experimental signatures of GRB in a kilometer-scale neutrino detector such as the proposed IceCube . We emphasize the importance of taking into account burst-to-burst fluctuations  as well as absorption of the neutrino signal in the Earth for both event rates and experimental signatures. Both effects produce additional and striking signatures with discriminating sensitivity to the value of the bulk Lorentz factor whose value is only indirectly inferred from other astronomical observations .

The observation of GRB neutrinos over a cosmological baseline has scientific potential beyond testing the “best-buy” fireball model: the observations can test with unmatched precision special relativity and the equivalence principle, and study oscillating neutrino flavors over the ultimate baseline of $`z1`$ .

## 2 Calculation of GRB Neutrino Rates and Signatures

In calculating the event rates and experimental signatures of GRB neutrinos in a high energy neutrino telescope we follow the model of Waxman and Bahcall  as implemented by Halzen and Hooper . We have normalized the neutrino flux to the energy rate injected in the Universe needed to explain the observed cosmic ray (CR) spectrum above $`10^{19}`$ eV, $`\dot{E}_{\mathrm{CR}}=4\times 10^{44}\mathrm{ergs}\mathrm{Mpc}^3\mathrm{yr}^1`$. This energy rate was calculated in reference , assuming a cosmological distribution of sources and taking into account CR propagation in the Cosmic Microwave Background Radiation (CMBR). The value quoted above corresponds to the “low redshift” ($`z<1`$) energy generation rate of CRs. Note that, because of the absorption of highest energy cosmic rays by CMBR photons, one could further increase $`\dot{E}_{\mathrm{CR}}`$ without directly affecting their observed flux. Waxman and Bahcall  have calculated an upper limit to the cosmic rate production rate in the whole Universe, assuming the energy generation rate evolves rapidly with redshift following the luminosity density evolution of quasi-stellar sources. They obtain an upper limit which is $`3\dot{E}_{\mathrm{CR}}`$. It is interesting however to mention that $`\dot{E}_{\mathrm{CR}}`$ is comparable to that produced in $`\gamma `$-rays by cosmological GRBs (which are not expected to be absorbed by the intervening backgrounds since the typical photon energy is below 1 MeV). Assuming the efficiency with which electrons (which ultimately produce the observed photons by synchrotron radiation) and protons is the same inside the GRB fireball, the value of $`\dot{E}_{\mathrm{CR}}`$ quoted above might well be closer to the actual value. For this, and other reasons, it is anyway unlikely that our calculations are accurate to better than a factor 3 or so. Moreover our neutrino event rate calculation might be conservative in this respect.

The neutrino flux is given by

$`{\displaystyle \frac{dN_\nu }{dE_\nu }}=\{\begin{array}{cc}\frac{A}{E_B}\frac{1}{E_\nu };E_\nu <E_B\hfill & \\ \frac{A}{E_\nu ^2};E_\nu >E_B\hfill & \end{array}`$ (1)
where $`A`$ is a normalization constant that is determined from energy considerations as explained above and in reference . The observed neutrino rates are calculated by folding this generic flux over the distributions of individual bursts in distance, $`f(z)`$, energy, $`f^{}(E_{\mathrm{GRB}})`$ and boost factor, $`f^{\prime \prime }(\mathrm{\Gamma })`$:

$$N_\nu \frac{dN_\nu }{dE_\nu }(E_{\mathrm{GRB}},z,\mathrm{\Gamma },E_\nu )P(E_\nu )f(z)f^{}(E_{\mathrm{GRB}})f^{\prime \prime }(\mathrm{\Gamma })𝑑E_{\mathrm{GRB}}𝑑z𝑑\mathrm{\Gamma }𝑑E_\nu .$$
(2)
where $`E_{\mathrm{GRB}}`$ is the energy emitted by a particular GRB, $`z`$ its redshift and $`\mathrm{\Gamma }`$ the boost factor. $`P(E_\nu )`$ represents the efficiency of detecting a neutrino of energy $`E_\nu `$. The first two distributions can be modelled after observations: a cosmological distribution in distance, and an energy distribution which assumes that ten percent of GRB produce more energy than average by a factor of ten, and one percent by a factor of 100 .

Most important however are the fluctuations in the $`\mathrm{\Gamma }`$ factor around its average phenomenological value of $`10^2`$$`10^3`$. The fluctuations in $`\mathrm{\Gamma }`$ affect the efficiency for producing pions in the $`p\gamma `$ collisions in the fireball as $`\mathrm{\Gamma }^4`$ , as well as the break energy, $`E_B`$, which varies as $`\mathrm{\Gamma }^2`$. Unfortunately the distribution in $`\mathrm{\Gamma }`$ cannot even be guessed at. Nevertheless, it is critical in making quantitative predictions . The physics is clear. In GRBs, high luminosities are emitted over short times, therefore the large photon density would render GRB opaque unless $`\mathrm{\Gamma }`$ is very large. Only transparent sources with large boost factors emit photons. They are however relatively weak neutrino sources because the actual photon target density in the fireball is diluted by the large Lorentz factor. An even moderately reduced value of $`\mathrm{\Gamma }`$ will produce a prolific neutrino source.

We remind the reader that the results obtained from Eq. 2 are at variance with the neutrino rate obtained by multiplying the average rate per burst by the number of bursts. Even neglecting the all important fluctuations in $`\mathrm{\Gamma }`$, there is no such concept as an average GRB. E.g. for $`\mathrm{\Gamma }=300`$, the correct computation of Eq. 2 yields a rate of $`75`$ events per km2 and year, roughly an order of magnitude larger than the prediction obtained by neglecting the observed burst-to-burst fluctuations in distance and energy. Another consequence of fluctuations is that the signal is dominated by a few very bright bursts, which greatly simplifies their detection.

As pointed out above, the average neutrino energy varies with the square of boost factor and therefore the calculated event rates, especially their dependence on $`\mathrm{\Gamma }`$, is strongly affected by the fact that higher energy neutrinos are preferentially absorbed in the Earth before reaching the detector . Also this effect has been neglected in all previous calculations. For instance, the 75 events just mentioned are reduced by a factor 3 by absorption. One should on the other hand remember that by oscillations, a large fraction of $`\nu _\mu `$ can oscillate into $`\nu _\tau `$ which penetrate the Earth . It is in this context important to realize that a kilometer-scale detector such as IceCube can measure the energy of the neutrinos. Therefore, signal events can be separated from the low energy atmospheric background by energy measurement, which allows the identifications of neutrinos from all directions and not just in the hemisphere where they pass through the Earth.

## 3 Results

We calculate the flux of neutrinos from GRB in the fireball model following reference . The number of protons in the fireball is fixed by the assumption that they are the source of the ultra high energy cosmic rays above $`10^{19}`$ eV. The results are shown in Figure 1 which shows the $`\nu _\mu +\overline{\nu }_\mu `$ flux from GRBs for different values of $`\mathrm{\Gamma }`$. The fluxes have been multiplied by $`E_\nu ^2`$ so that they represent the energy emitted in the form of neutrinos. Notice the variation of the break in the spectrum $`E_B`$ as $`\mathrm{\Gamma }^2`$. Also, for values of $`\mathrm{\Gamma }`$ below $``$ 100 the fireball becomes opaque to protons and at this point the total amount of energy available for neutrino production is converted. The neutrino flux, which roughly scales as $`\mathrm{\Gamma }^4`$, saturates and no longer grows with decreasing boost factor; see Figure 1.

As usual , we calculate the event rate per GRB by folding the neutrino flux with the probability of detecting a muon produced in a $`\nu _\mu `$ or $`\overline{\nu }_\mu `$ interaction. Fluctuations in distance $`z`$ and total GRB energy, or fluency, $`E_{\mathrm{GRB}}`$ are accounted for by a Monte Carlo simulation we have developed for this purpose. We simulate a large number of GRBs at different zenith angles assuming an isotropic distribution and subsequently obtain the event rate per year by assuming 1000 GRBs/year. Absorption of the neutrinos in the Earth prior to reaching the detector is taken into account using the density profile of reference . It is also important to implement the fact that above the detector there is a limited column density of atmosphere and ice available for neutrino detection.

In Table 1 we separately show, for different $`\mathrm{\Gamma }`$’s, the event rate of upgoing, i.e. neutrinos that cross the Earth before interacting near the detector, and downgoing neutrinos. The first two columns show for comparison the number of events when neither absorption nor the limited amount of target above the detector are taken into account. It is clear that both effects play an important role in obtaining the correct event rate. This is not surprising: the Earth becomes opaque to neutrinos of energy around 100 TeV, and the muon range exceeds the $``$ 2 km vertical depth of the IceCube detector at energies around $`1`$ TeV.

The dependence of the number of events on the $`\mathrm{\Gamma }`$ factor is shown in Figure 2. Two competing effects determine the shape of the curve. The event rate decreases with increasing $`\mathrm{\Gamma }`$ following the dependence of the neutrino flux which varies as $`\mathrm{\Gamma }^4`$. This decrease is partially offset because higher energy neutrinos resulting from larger boost factors are more efficiently detected. For low values of $`\mathrm{\Gamma }`$, below about 100, the saturation of the total energy available for neutrino production is seen. On the other side of the $`\mathrm{\Gamma }`$ range, for large values of $`\mathrm{\Gamma }`$, the spectrum is very flat ($`E^1`$) up to $`70`$ PeV where the absorption by the Earth dominates. This reduces the event rate of upgoing events for values of $`\mathrm{\Gamma }`$ above 1000 as can be seen in Figure 2. In the end the downgoing and upgoing event rates are similar except for large values of $`\mathrm{\Gamma }`$. It is important to keep in mind that although downgoing neutrinos are not affected by absorption, their detection is limited by the column density of matter available for neutrino interaction.

In Figure 3, we show the energy dependence of upgoing and downgoing $`\nu _\mu +\overline{\nu }_\mu `$ events for three representative values of $`\mathrm{\Gamma }`$. The zenith angle distribution of upgoing neutrinos is shown in Figure 4, i.e. $`1<\mathrm{cos}(\theta _{\mathrm{zenith}})<0`$. As $`\mathrm{\Gamma }`$ increases, the higher energy neutrinos are attenuated by the Earth. This explains why, as $`\mathrm{\Gamma }`$ increases, the distributions increasingly resemble an exponential attenuation function.

### 3.1 $`\nu _\tau +\overline{\nu }_\tau `$ events

$`\nu _\tau `$ production is expected to be very small in GRBs and in general in any astrophysical environment where $`\nu `$’s are produced in $`pp`$ or $`p\gamma `$ collisions. Several calculations suggest a ratio :

$$F_{\nu _\tau /\nu _\mu }=\frac{\mathrm{\Phi }(\nu _\tau +\overline{\nu }_\tau )}{\mathrm{\Phi }(\nu _\mu +\overline{\nu }_\mu )}10^5,$$
(3)
Oscillation scenarios in which $`\nu _\mu `$’s convert into $`\nu _\tau `$’s can however provide abundant sources of $`\nu _\tau `$’s. Assuming a typical value $`F_{\nu _\tau /\nu _\mu }=0.5`$ suggested by SuperKamiokande measurements , we obtain the double bang  111Events in which two separated showers can be identified, one initiated by the struck nucleon and the other by the decay of the $`\tau `$ produced in the $`\nu _\tau `$ charged current interaction. event rates shown in Table 2. The probability of detecting a $`\nu _\tau `$ induced double bang event in IceCube is typically two orders of magnitude smaller than the probability of observing a $`\nu _\mu `$ at energies around 10 PeV . This explains the small values of the event rates in Table 2. In the event rate calculation we have accounted for the energy loss of the $`\nu _\tau `$’s when they propagate along the Earth’s interior. This produces a pileup of events around 100 TeV, as pointed out in , reducing the number of upgoing double bang events with respect to the downgoing ones. This is due to the probability of $`\nu _\tau `$-induced double bang detection which is limited to a broad peak between $``$ 1 PeV and $``$ 100 PeV (outside of which it is negligible). This also explains the fact that the energy distribution of the event rate peaks in the vicinity of 10 PeV (see Figure 5).

We also calculated the $`\nu _\tau +\overline{\nu }_\tau `$ events that would be detected by the appearance of a $`\tau `$ which decays to $`\mu `$ just below the detector. The event rates are also shown in Table 2. In this case the energy distribution of the events peaks around the energy at which the events pileup (100 TeV); see Figure 5. The probability of detecting the $`\mu `$ is $`17\%`$ of the probability of detecting it in a $`\nu _\mu `$ interaction due to the branching ratio of the $`\tau `$ to $`\mu `$ decay channel. This accounts for a factor $``$ 6 difference between $`\nu _\tau \tau \mu `$ in Table 2 and the $`\nu _\mu +\overline{\nu }_\mu `$ rate in Table 1.

The expected event rates generated by both mechanisms (double bang and $`\tau \mu `$ decay) have very different and characteristic zenith angle distributions. These are shown in Figure 6. For large column depths inside the Earth (i.e. $`\mathrm{cos}(\theta _{\mathrm{zenith}})1`$), the events pileup around 100 TeV and the small probability of detecting a double bang event at this energy reduces the number of double bangs with respect to events in which a $`\mu `$ is detected. When $`\mathrm{cos}(\theta _{\mathrm{zenith}})0`$, because the amount of matter neutrinos have to cross is very small, there is no pileup of events and the angular distributions have roughly the same shape. Figure 6 also shows the zenith angle distribution of events produced by $`\nu _\mu +\overline{\nu }_\mu `$ for comparison, making it clear that despite of the flatter distribution of the $`\mu `$’s from $`\nu _\tau +\overline{\nu }_\tau `$ they are still outnumbered by $`\mu `$’s from $`\nu _\mu `$ interactions.

### 3.2 Fluctuations in $`\mathrm{\Gamma }`$

Finally, we discuss burst to burst fluctuations in the boost factor $`\mathrm{\Gamma }`$. So far we, quite unrealistically, assumed that $`\mathrm{\Gamma }`$ is single valued. We have, at present, no information on the distribution of $`\mathrm{\Gamma }`$ factors. It has been shown  that the neutrino rates can be significantly enhanced by fluctuations around the average value. It may be, of course, that the fluctuations in energy and distance which we took into account following the experimental evidence, already reflect some or all of the fluctuations in the boost factor. Following reference  we will illustrate the effect by assuming Gaussian distributions with half width $`\sigma `$. The results are shown in Table 3 for different values of $`\mathrm{\Gamma }`$. It is interesting to note that for $`\mathrm{\Gamma }=100`$ the event rates are almost independent of the value of $`\sigma `$. For $`\mathrm{\Gamma }`$ around 100 the event rate is weakly dependent on $`\mathrm{\Gamma }`$ due to the saturation of the amount of energy that goes into neutrino production (see Figures 1 and 2) and hence fluctuations in the value of the boost factor do not affect the event rate. This is not the case when $`\mathrm{\Gamma }=300`$ or 1000 since the event rate behaves roughly as $`\mathrm{\Gamma }^4`$ in that $`\mathrm{\Gamma }`$ range. Figure 7 shows more clearly the dependence of the event rate on $`\sigma `$ for different values of $`\mathrm{\Gamma }`$ illustrating these points. It can be easily shown that when $`\sigma \mathrm{\Gamma }`$ the following relation holds:

$$\frac{\mathrm{Rate}(\mathrm{\Gamma }\pm \sigma )}{\mathrm{Rate}(\mathrm{\Gamma })}\frac{\sigma }{\mathrm{\Gamma }},$$
(4)
explaining why for $`\mathrm{\Gamma }=1000`$ the variation of the event rate with $`\sigma `$ is not as strong as when $`\mathrm{\Gamma }=300`$, even though the rate scales with $`\mathrm{\Gamma }^4`$ in both cases. The same comments apply to double bang events produced by $`\nu _\tau `$’s and events in which the produced $`\tau `$ decays to $`\mu `$ (both are shown in Table 4).

## 4 Summary and conclusions

We have investigated potential signatures in $`\mathrm{km}^3`$ telescopes of high energy neutrino fluxes produced in $`p\gamma `$ interactions in GRB environments. We stress the fact that the rate is dominated by fluctuations in distance, GRB energy and in the bulk Lorentz factor $`\mathrm{\Gamma }`$ of the expanding GRB fireball. We have used an euclidean distribution to account for fluctuations in distance. Using a distribution in which GRB’s follow the star formation rate, as suggested in , the event rates are reduced only by $`20\%`$. On the other hand, a cosmological distribution following galaxies  allows for more close GRB’s so that the event rate increases by a factor of 20 (see Table 5). Neutrino telescopes may help to constrain the distribution in distance of GRBs. The most relevant parameter is $`\mathrm{\Gamma }`$, which determines the rate of $`p\gamma `$ interactions and hence the amount of energy that goes into neutrino production. For small values of $`\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }50`$ or less) the expansion of the GRB fireball is not sufficiently fast and the large photon density makes it opaque to $`p\gamma `$, efficiently producing pions — the parents of the neutrinos — and saturating the amount of energy available for neutrino production. Due to the scaling of the energy break in the spectrum with $`\mathrm{\Gamma }^2`$, mostly low energy neutrinos are produced whose detection efficiency is smaller due to the small range of the muon at low energies and the relatively high energy threshold of the neutrino telescopes ($`100`$ GeV). However the large neutrino flux compensates for the small detection efficiency. For large values of $`\mathrm{\Gamma }`$ ($`\mathrm{\Gamma }>300`$), an increasingly larger fraction of the neutrino energy goes into the high energy part of the spectrum, however the overall amount of energy is very small, producing small event rates. Values of $`\mathrm{\Gamma }<100`$ give rise to large number of events that would even be observable in existing smaller neutrino telescopes such as AMANDA  (one has to scale the results down by roughly two orders of magnitude to account for the smaller effective area of the detector).

We have shown that absorption of the upgoing $`\nu _\mu `$’s inside the Earth, as well as the limited column of matter available for downgoing neutrino interactions, play a relevant role, making upgoing and downgoing event rates roughly equal. The change of the energy break of the spectrum with $`\mathrm{\Gamma }^2`$ combined with the absorption of the Earth is reflected in the zenith angle distributions of the event rates which may give some complementary information about $`\mathrm{\Gamma }`$.

GRB neutrino detection with $`\mathrm{km}^3`$ neutrino telescopes also has the potential to investigate $`\nu _\mu \nu _\tau `$ oscillations over a baseline of 1000 Mpc. Double bang $`\nu _\tau `$ induced events offer an unmistakeable signature which allows downgoing $`\nu _\tau `$ detection. A $`\mathrm{km}^3`$ telescope operating for 10 years may detect $``$ 10 downgoing double bang events if $`\mathrm{\Gamma }=100`$ without any potential background. Upgoing double bang events are not going to be detected since the $`\nu _\tau `$’s pileup around 100 TeV, where the probability of a double bang detection is negligible. Muons produced in $`\tau `$ decays are outnumbered by $`\mu `$’s from $`\nu _\mu `$ interactions at least for the type of fluxes expected from GRBs ; besides, their energy distribution does not show a clear and characteristic signature so they are difficult to identify.

In summary, neutrino telescopes open up the possibility of determining the value of $`\mathrm{\Gamma }`$ and its fluctuations, as well as the possibility of identifying $`\nu _\tau `$’s. They are powerful instruments to reveal important astrophysical information about the most energetic objects ever observed in the universe and about neutrino oscillation scenarios over cosmological baselines.

## Acknowledgments

This research was supported in part by the US Department of Energy under grant DE-FG02-95ER40896 and in part by the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation. J.A. thanks the Department of Physics, University of Wisconsin, Madison and the Fundación Caixa Galicia for financial support.