text:
# Enhanced magnetocaloric effect in frustrated magnets

## I Introduction

The magnetocaloric effect consists in heating or cooling of matter in response to variations of external magnetic field. The technique of adiabatic demagnetization based on the magnetocaloric effect has been successfully utilized over the years to reach temperatures in a sub-Kelvin range. Lounasmaa  Because of its technological simplicity and independence of gravity the magnetic cooling steadily attracts attention as a rival to the dilution refrigerators, especially in satellite applications. Hagmann ; McRae ; Kushino  Recent experimental studies of the magnetocaloric effect in specially designed materials have revealed the potential applicability of the magnetic cooling technique for room temperature refrigeration as well. Pecharsky ; Giguere ; Tegus

Paramagnetic salts, which are standard refrigerant materials for the low-temperature magnetic cooling, contain dilute systems of magnetic dipoles. The entropy of an ideal paramagnet depends on magnetic field $`H`$ and temperature $`T`$ via $`H/T`$. Therefore, during adiabatic demagnetization, temperature of a paramagnet decreases linearly with a field: $`T/H=\mathrm{const}`$. The higher the density of the magnetic moments and their spin number, the greater the cooling power of a refrigerant is. With increased density of spins, however, the strength of residual interactions between magnetic dipoles grows and leads to an ordering or a spin glass transition. The transition temperature $`T_c`$, which is roughly given by a strength of magnetic interactions, limits the lowest temperatures achievable with paramagnetic salts. Even above $`T_c`$ the cooling becomes less efficient due to a reduced entropy, and the adiabatic temperature change is slower than in an ideal paramagnet. The cooling rate $`(T/H)_S`$ of a salt does not, therefore, exceed the corresponding rate of a paramagnet: $`(T/H)_S^{\mathrm{para}}=T/H`$.

The experimental and theoretical studies over the past decade have established a new distinct class of magnetic materials called geometrically frustrated magnets.Ramirez ; Greedan ; Moessner  Magnetic ions in these systems form special types of lattices, Fig. 1. Despite interaction between neighboring spins, strongly frustrated magnets remain disordered and possess finite entropy at temperatures well below the Curie-Weiss constant. This property suggests frustrated magnets as prospective candidates for use in the adiabatic demagnetization refrigerators. In fact, gadolinium gallium garnet Gd3Ga5O12 (GGG) has long been known as a suitable refrigerant material. GGGcool  Recently, GGG with its unique phase diagram in magnetic field attracted a lot of interest from the point of view of geometric magnetic frustration on a garnet lattice.GGGphd  In the present work we demonstrate a close relation between an enhanced magnetocaloric effect and strong geometric frustration. In particular, $`(T/H)_S`$ in a frustrated magnet can exceed the cooling rate of an ideal paramagnet by more than an order of magnitude.

Specifically, we study an unusual magnetothermodynamics of classical Heisenberg antiferromagnets on several frustrated lattices. The choice of classical models is motivated by large values of spins of the magnetic ions in frustrated magnets: Cr3+ ($`S=3/2`$), Fe3+ ($`S=5/2`$), and Gd3+ ($`S=7/2`$).Ramirez ; Greedan  The studied frustrated structures include two- (2D) and three-dimensional (3D) networks of corner-sharing triangles known as a kagome, Fig. 1a, and a garnet lattice, Fig. 1b, and a 3D network of corner-sharing tetrahedra of a pyrochlore lattice, Fig. 1c. In all three cases the magnetic Hamiltonian includes exchange interaction of strength $`J`$ between the nearest-neighbor classical spins of unit length and the Zeeman energy:

$$\widehat{}=J\underset{ij}{}𝐒_i𝐒_j𝐇\underset{i}{}𝐒_i.$$
(1)
The main distinction of strongly frustrated magnets is a large, macroscopic degeneracy of their classical ground states in zero magnetic field. It arises due to a so called under-constraint: condition of the minimal energy for every elementary frustrated unit, triangle or tetrahedron, is not sufficient to fix all microscopic degrees of freedom. Ramirez ; Moessner  In applied magnetic field the macroscopic degeneracy persists up to the saturation field $`H_{\mathrm{sat}}`$, which is $`H_{\mathrm{sat}}=6J`$ for kagome and garnet antiferromagnets and $`H_{\mathrm{sat}}=8J`$ for a pyrochlore antiferromagnet. fafm ; kagome  At zero temperature frustrated magnets loose their degeneracy and become completely polarized above $`H_{\mathrm{sat}}`$. Transformation from the unique ground state at $`H>H_{\mathrm{sat}}`$ to an infinitely degenerate ground state at $`H<H_{\mathrm{sat}}`$ is accompanied by condensation of a macroscopic number of zero-energy modes. kagome ; Honecker  This should be contrasted with the behavior of non- or weakly-frustrated magnets, where the antiferromagnetic state below $`H_{\mathrm{sat}}`$ is described by a certain wave-vector and, therefore, the phase transition at $`H=H_{\mathrm{sat}}`$ corresponds to condensation of one or a few zero-energy modes. In the vicinity of $`H_{\mathrm{sat}}`$ the thermodynamic properties of a strongly frustrated magnet are, therefore, controlled by flat branches of excitations with vanishing energy at $`H=H_{\mathrm{sat}}`$.

The following theoretical analysis is divided into two parts. The next section gives a general thermodynamic description of magnets with a macroscopic number of soft excitations. We also calculate the number of such modes for the three frustrated spin models. The subsequent section describes Monte Carlo results for various thermodynamic quantities of frustrated models and also simulations of the adiabatic demagnetization process for a few frustrated and non-frustrated magnets. These are followed by discussion of suitable magnetic materials, which can exhibit the predicted behavior.

## II Theory

Our aim in this section is to understand the effect of soft modes on a cooling rate of a strongly frustrated magnet under adiabatic demagnetization. The cooling rate is related to isothermal characteristics by

$$\left(\frac{T}{H}\right)_S=T\frac{(S/H)_T}{C},$$
(2)
where $`C`$ is the specific heat in a constant magnetic field. We shall also use the normalized rate dividing (2) by the cooling rate of an ideal paramagnet $`(T/H)_S^{\mathrm{para}}=T/H`$. Values $`(T/H)_S^{\mathrm{norm}}>1`$ indicate that a magnet cools down faster than a paramagnet in a certain range of fields and temperatures. The cooling rate is enhanced for systems with large negative values of $`(S/H)_T`$.

### II.1 Magnetothermodynamics

A classical magnet on an $`N`$-site lattice has $`2N`$ degrees of freedom in total with two modes per spin. In the saturated phase above $`H_{\mathrm{sat}}`$ all $`2N`$ magnetic excitations are gapped. In the harmonic approximation the excitation energies of non-frustrated magnets are given by

$$\omega (𝐤)=HH_{\mathrm{sat}}+\epsilon (𝐤),$$
(3)
where $`\epsilon (𝐤)`$ is nonnegative and vanishes at one or a few points in the Brillouin zone. Non-Bravais lattices have several modes for every wave-vector $`𝐤`$ and below we imply summation over all such branches. At low temperatures anharmonic effects are small, which allows to find explicitly the fluctuation contributions to various physical quantities. In particular, the specific heat per one spin is $`C=1`$ ($`k_B1`$) and

$$\left(\frac{S}{H}\right)_T=\frac{1}{2}\underset{𝐤}{}\frac{1}{HH_{\mathrm{sat}}+\epsilon (𝐤)}.$$
(4)
Since $`\epsilon (𝐤)k^2`$ near its zeroes, $`(S/H)_T`$ at $`H=H_{\mathrm{sat}}`$ is a $`T`$-independent constant in 3D and has a weak logarithmic singularity in 2D.

The excitation spectra of strongly frustrated magnets in the saturated phase are special in that $`\epsilon (𝐤)0`$ in the whole Brillouin zone for several flat branches of excitations. We denote the number of such modes as $`N_4`$, $`N_2`$ being the number of ordinary dispersive modes (3) and $`N_2+N_4=2N`$. As field approaches $`H_{\mathrm{sat}}`$, the contribution from soft dispersionless excitations grows and the harmonic approximation fails. To study the effect of soft modes one needs to consider an anharmonic energy functional of the type:

$$E=\underset{𝐤}{}(HH_{\mathrm{sat}})|y_𝐤|^2+\underset{𝐤,𝐩,𝐪}{}V_{\mathrm{𝐤𝐩}}^𝐪y_{𝐤+𝐪}^{}y_{𝐩𝐪}^{}y_𝐩y_𝐤,$$
(5)
where $`y_𝐤`$’s denote collective coordinates of soft modes and $`V_{\mathrm{𝐤𝐩}}^𝐪J`$ are nonnegative interaction constants. Cubic terms vanish in the expansion by a spin-rotational symmetry. Also, an interaction between soft and the harmonic modes can be neglected at low temperatures. The above energy functional determines scaling of the corresponding contribution to the free energy:

$$F_4=\frac{N_4}{4}T\mathrm{ln}TTf\left(\frac{HH_{\mathrm{sat}}}{\sqrt{T}}\right).$$
(6)
The function $`f(u)`$ has the following asymptotes: $`f(u)f(0)+uf^{}(0)`$ with $`f(0)>0`$ and $`f^{}(0)<0`$ as $`u0`$, and $`f(u)(N_4/2)\mathrm{ln}u`$ as $`u\mathrm{}`$.

It follows from Eq. (6) that the specific heat and the magnetic susceptibility depend on magnetic field and temperature only via combination $`(HH_{\mathrm{sat}})/\sqrt{T}`$.kagome  In particular, the specific heat per one spin reaches a universal fractional value at $`H=H_{\mathrm{sat}}`$: $`C=(\frac{1}{2}N_2+\frac{1}{4}N_4)/N=1N_4/(4N)`$. For the magnetocaloric effect, Eq. (6) gives at $`H=H_{\mathrm{sat}}`$:

$$\left(\frac{S}{H}\right)_T\frac{N_4}{\sqrt{T}},$$
(7)
which significantly exceeds the contribution of dispersive modes Eq. (4). Thus, condensation of a macroscopic number of zero-energy modes at the saturation field greatly enhances the magnetocaloric effect in strongly frustrated magnets at low temperatures. An important characteristics of a magnet is, then, the number of soft modes $`N_4`$.

### II.2 Soft modes

#### II.2.1 Kagome antiferromagnet

In order to find the number of soft modes we perform classical spin-wave calculations in the saturated collinear state. A primitive unit cell of a kagome lattice contains one triangular plaquette. The corresponding Bravais lattice has $`N/3`$ states. Expanding spin components along the field in small transverse deviations we obtain a classical spin-wave Hamiltonian to a desired order. In the harmonic approximation two transverse polarizations for the spin excitations are independent. As a result, the spectrum consists of three double-degenerate modes:

$`\omega _1(𝐤)`$ $`=`$ $`H6J,`$
$`\omega _{2,3}(𝐤)`$ $`=`$ $`H3J\pm J\sqrt{3(1+2\gamma _𝐤)},`$ (8)
where $`\gamma _𝐤=\frac{1}{3}(\mathrm{cos}k_x+2\mathrm{cos}\frac{1}{2}k_x\mathrm{cos}\frac{\sqrt{3}}{2}k_y)`$. The total number of soft quartic modes at $`H_{\mathrm{sat}}=6J`$ is, therefore, $`N_4^{\mathrm{kag}}=2N/3`$. This is exactly twice the number of hexagon voids on a kagome lattice. The correspondence between such voids and soft modes has first been found for a kagome antiferromagnet in zero field Chalker  and was, later, extended to finite magnetic fields.kagome  Local soft modes in the collinear state at $`H=H_{\mathrm{sat}}`$ consist in alternate tilting of spins around elementary hexagons in two different directions perpendicular to the field. The geometric picture of local soft modes can be straightforwardly generalized for a classical spin model on a checker-board lattice, which is the other planar frustrated lattice with voids around empty squares.Honecker  The checker-board antiferromagnet has $`N`$ local zero-energy modes at $`H_{\mathrm{sat}}=8J`$.

#### II.2.2 Pyrochlore antiferromagnet

The spin-wave calculations for a pyrochlore antiferromagnet yield four double-degenerate modes defined in the Brillouin zone of a face-centered cubic (fcc) lattice, which has $`N/4`$ states:

$`\omega _{1,2}(𝐤)`$ $`=`$ $`H8J,`$
$`\omega _{3,4}(𝐤)`$ $`=`$ $`H4J\pm 2J\sqrt{1+3\nu _𝐤}`$ (9)
with $`\nu _𝐤=\frac{1}{6}(\mathrm{cos}k_x+\mathrm{cos}k_y+4\mathrm{cos}k_z\mathrm{cos}\frac{1}{2}k_x\mathrm{cos}\frac{1}{2}k_y)`$. The total number of zero modes in a pyrochlore antiferromagnet at the saturation field is $`N_4^{\mathrm{pyr}}=N`$, i.e. the same as for the checker-board model. The geometric interpretation of these modes becomes, however, more complicated. A dispersionless soft mode can still be attributed to a local perturbation in the collinear structure: alternate tilting of spins around a closed line, which shares two vertices with each crossed tetrahedron. The shortest such loop on a pyrochlore lattice is a hexagon. Simple counting yields $`N`$ hexagons for an $`N`$-site pyrochlore lattice. If a one-to-one correspondence between the soft modes and the hexagons holds, then the above number would imply $`2N`$ soft modes, which is twice more than the spin-wave prediction. The contradiction is resolved in favor of the spin-wave result by establishing that some of the hexagon loops on a pyrochlore lattice are linearly dependent with total $`N/2`$ linear relations. This brings the number of independent hexagon spin-loops (modes) to $`N/2`$ in agreement with the spin-wave expansion. Below, we briefly sketch a proof of the above statement.

A pyrochlore lattice, Fig. 1c, is formed with two types of tetrahedra. To simplify geometric representation we draw a pyrochlore structure as an fcc lattice constructed from one type of tetrahedra. A cut through such a lattice along $`(001)`$ plane is shown in Fig. 2. Hexagons in Fig. 2 correspond to closed loops which include two tetrahedra in one layer and one tetrahedron in an adjacent layer. Soft spin modes on two such loops are shown in Fig. 2 with pluses and minuses, which indicate alternate tilting of spins from the field direction. The first type of a linear relation between spin-loops appears for a block formed by four tetrahedra numbered 1–4 in Fig. 2a. A superposition of the local soft modes corresponding to hexagon loops $`l_{123}`$ and $`l_{134}`$ leaves two inner spins unmoved and tilts only outer spins around the block perimeter. The same eight-spin mode is obtained by a superposition of spin-loops $`l_{124}`$ and $`l_{234}`$. Thus, every such block of four tetrahedra provides one linear relation between four loops. This reduces the number of independent spin-loops by the number of such blocks $`N/4`$. The second type of a linear relation corresponds to a block of six tetrahedra: four in one layer, 1, 3, 5, and 6, and two in adjacent layers, 2 and 7, Figs. 2a,b. A superposition of spin-loops $`l_{123}`$ and $`l_{256}`$ corresponds to the same mode as a superposition of $`l_{157}`$ and $`l_{367}`$. Counting all such blocks brings additional $`N/4`$ linear relations. This makes in total $`N/2`$ linear relations between hexagonal spin-loops and, hence, $`N_4^{\mathrm{pyr}}=2(NN/2)=N`$ soft modes for a pyrochlore antiferromagnet at the saturation field. Composite low-energy excitations on six-spin loops, which are similar to the above soft modes, have been recently observed in neutron scattering studies of pyrochlore antiferromagnet ZnCr2O4 in zero field.Broholm

#### II.2.3 Garnet antiferromagnet

The structure of a garnet lattice is by far the most complicated among frustrated lattices. In GGG the magnetic Gd3+ ions form two interpenetrating triangular sublattices with 24 atoms per a cubic unit cell. Exchange interaction couples each spin with its nearest-neighbors from the same sublattice. Spins on two different sublattices interact only via weaker dipolar forces. In the present study we focus on the behavior of an exchange model and, therefore, consider only one sublattice of corner-sharing triangles with 12 spins per a unit cell represented, e.g., by dark circles in Fig. 1b. The nearest-neighbor exchange Hamiltonian on a garnet lattice can be written as

$``$ $`=`$ $`J{\displaystyle \underset{i}{}}[𝐒_{1i}(𝐒_{2i}+𝐒_{3i}+𝐒_{4i+x}+𝐒_{5i+x})+𝐒_{2i}(𝐒_{3i}+𝐒_{6i+y}+𝐒_{7i+y})+𝐒_{3i}(𝐒_{8i+z}+𝐒_{9i+z})`$ (10)
$`+𝐒_{4i}(𝐒_{5i}+𝐒_{8i+z}+𝐒_{10i+z})+𝐒_{5i}(𝐒_{6ix}+𝐒_{11i})+𝐒_{6i}(𝐒_{7i}+𝐒_{11i+x})+𝐒_{7i}(𝐒_{9iy}+𝐒_{12i})`$
$`+𝐒_{8i}(𝐒_{9i}+𝐒_{10i})+𝐒_{9i}𝐒_{12i+y}+𝐒_{10i}(𝐒_{11i}+𝐒_{12i})+𝐒_{11i}𝐒_{12i}]+_{\mathrm{Zeeman}},`$
where the first spin index $`j=1`$–12 is the site number in a unit cell and the second spin index $`i`$ is the cell number on a cubic lattice; $`i\pm x`$,… being adjacent cells along three orthogonal directions.

After Fourier transformation the spin-wave energies are given by eigenvalues of a $`12\times 12`$ matrix. The corresponding equation factorizes into a product of $`(\omega H+6J)^4`$ and an eighth-order polynomial in $`\omega `$. math  At $`H=6J`$ there are exactly four zero roots for an arbitrary value of the wave-vector. Zero eigenvalues determine $`N_4^{\mathrm{gar}}=2N/3`$ soft modes at the saturation field, which is the same number as for a kagome antiferromagnet. Thus, the term hyper-kagome used sometimes to refer to a garnet structure receives an extra justification. Geometric picture of soft modes in a hyper-kagome antiferromagnet is, however, much more complicated than in its planar counterpart. The shortest closed loop on a garnet lattice, which corresponds to an elementary soft mode, crosses ten triangles. At the same time each bond between a pair of nearest-neighbor spins participates in five such loops. private  This should be contrasted with simple hexagon loops on a kagome lattice, where every bond is included in one loop only. The total number of ten-triangle loops for a garnet lattice is calculated to be $`N`$. A comparison with $`N/3`$ soft modes for one transverse polarization found in the spin-wave analysis suggests that there are $`2N/3`$ linear relations between closed loops for a garnet structure. An underlying geometric picture behind these relations remains to be clarified.

To conclude this section we note that a pyrochlore lattice is the most frustrated one among the above three structures, because it has the largest number of soft modes. The same is also true in zero magnetic field, where a classical Heisenberg pyrochlore lattice antiferromagnet remains disordered down to lowest temperatures, Moessner\_Chalker  whereas a classical kagome antiferromagnet develops a unique planar triatic order parameter. kagome ; Chalker  As a result, a pyrochlore antiferromagnet should exhibit the fastest cooling rate under adiabatic demagnetization among frustrated magnets.

## III Monte Carlo simulations

Monte Carlo (MC) simulations of the three spin models have been performed with the standard Metropolis algorithm. To improve statistics at low temperatures we introduce a maximum change $`\mathrm{\Delta }S^z=3T`$ for a choice of a new spin direction in the local coordinate frame. In this way an acceptance rate for spin moves is kept at the level of 30%. (Other values $`\mathrm{\Delta }S^z/T1`$–5 will work as well.) MC simulations are done on finite clusters with periodic boundary conditions. Typical cluster sizes are $`N=3\times 18^2=972`$ for a kagome lattice, $`N=4\times 6^3=864`$ for a pyrochlore lattice, and $`N=12\times 5^3=1500`$ for a garnet lattice. Comparison with smaller and larger clusters have been made to insure that finite-size corrections are small. Typical MC runs are carried with $`10^5`$ MC steps per spin for equilibration followed by averaging over $`10^5`$ measurements made in intervals of 5–10 MC steps. Statistical errors are estimated by performing several ($`5`$) runs from different initial configurations. The errors for all presented data are small and do not exceed symbol sizes. The specific heat is calculated from the variance of the total energy: $`C=(E^2E^2)/T^2`$, while the magnetocaloric effect is expressed via a cumulant of the energy and the magnetization: $`(S/H)_T=(EMEM)/T`$.

We present in Fig. 3 the specific heat per one site of kagome and pyrochlore antiferromagnets plotted with respect to the scaling parameter $`u=(HH_{\mathrm{sat}})/\sqrt{T}>0`$. The data have been collected during field-scans at fixed temperatures. Points from different scans, except of the highest temperature of $`T=0.1J`$, fall on the same curve confirming the predicted scaling form. The specific heat at $`H=H_{\mathrm{sat}}`$ ($`u=0`$) reaches values which are close to the asymptotic results derived from the number of soft quartic modes: $`C=5/6`$ and $`C=3/4`$ for a kagome and a pyrochlore antiferromagnet, respectively. The magnetocaloric effect near saturation is presented in Fig. 4. The main plot shows temperature dependence of $`(S/H)_T`$ exactly at the saturation field for kagome, pyrochlore and square lattice (576-site cluster) antiferromagnets. Numerical results for a garnet antiferromagnet are not included because at $`H>H_{\mathrm{sat}}`$ they virtually coincide with the data for a kagome antiferromagnet. Significant difference (Fig. 4) between the frustrated magnets on one side and a non-frustrated square lattice antiferromagnet on the other side confirms an enhancement of the magnetocaloric effect related to magnetic frustration. A range of magnetic fields in the vicinity of $`H_{\mathrm{sat}}`$, where such an enhancement occurs, is illustrated on the inset.

Direct comparison of the cooling rates of magnets with and without frustration is made in Fig. 5, where results for a pyrochlore and a square lattice antiferromagnet are plotted for several field scans. The two magnets have the same saturation field $`H_{\mathrm{sat}}=8J`$. For a pyrochlore antiferromagnet this field corresponds to a crossover from the nondegenerate saturated phase into a highly degenerate state below $`H_{\mathrm{sat}}`$. The crossover is accompanied by a large increase in the adiabatic temperature variation, which exceeds the cooling rate of an ideal paramagnet by more than an order of magnitude. In contrast the spin model on a square lattice has a real transition near $`H_{\mathrm{sat}}`$ into a nondegenerate antiferromagnetic state with divergent correlation length. The cooling rate becomes negative below the transition, which means that a magnet heats up during demagnetization rather than cooling down. There is also no significant enhancement of the cooling rate above $`H_{\mathrm{sat}}`$ apart from a small narrow peak due to an anomaly in $`(S/H)_T`$ at the transition.

Finally, a full adiabatic demagnetization process can be simulated with MC by measuring at fixed $`T`$ and $`H`$ isothermal characteristics $`C`$ and $`(S/H)_T`$ to find the cooling rate (2) and, then, choosing a new value of temperature according to $`\mathrm{\Delta }T_{ad}/\mathrm{\Delta }H=(T/H)_S`$ for sufficiently small $`\mathrm{\Delta }H`$. Adiabatic temperature changes for demagnetization from two starting temperatures $`T=J`$ and $`T=0.1J`$ and the same initial field $`H=9J`$ are shown for garnet, pyrochlore, and cubic ($`N=512`$) antiferromagnets in Fig. 6. Since saturation fields for the three magnets are different, we use rescaled exchange constants $`J^{}=0.75J`$ for a pyrochlore and $`J^{}=0.5J`$ for a cubic antiferromagnet to set the common saturation field of $`H_{\mathrm{sat}}=6J`$ in all three cases. During demagnetization the non-frustrated cubic antiferromagnet has the smallest relative temperature variation and starts to heat up below the transition into antiferromagnetically ordered state. The adiabatic temperature decrease for frustrated magnets is much larger with the largest effect exhibited by a pyrochlore antiferromagnet. Note, that temperature of a pyrochlore antiferromagnet drops by more than ten times as field decreases from $`H_{\mathrm{ini}}=1.5H_{\mathrm{sat}}`$ to $`H_{\mathrm{fin}}=0.8H_{\mathrm{sat}}`$. To achieve a similar cooling with a paramagnetic salt magnetic field has to be decreased at least by a factor of ten. The two distinctive features of the demagnetization of a pyrochlore antiferromagnet are (i) a persistent temperature decrease up to a zero field, though with a reduced cooling rate below $`H_{\mathrm{sat}}`$, and (ii) an increase in the relative temperature variation as a starting temperature goes down, which is a consequence of the low-temperature singularity (7).

## IV Discussion

The practical limitation on experimental observation of an enhanced magnetocaloric effect comes from a strength of the saturation field. At present, only frustrated magnets with exchange constants of the order of few Kelvin or less can be saturated with dc magnetic fields. The prime example is gadolinium gallium garnet with $`H_{\mathrm{sat}}1.7`$ T,GGGphd  which is a popular refrigerant material in an applied experimental research on the magnetic cooling. GGGcool ; Kushino  Quantitative comparison of the cooling rates of an ideal paramagnet and a Heisenberg garnet antiferromagnet is made in Fig. 6. Starting with the same temperature a garnet antiferromagnet cools down to much lower temperatures than a paramagnet in a wide range of magnetic fields. There are also two pyrochlore antiferromagnets Gd2Ti2O7 and Gd2Sn2O7,Raju ; Ramirez02 ; Bramwell  which have saturation fields about 7 T. It is an interesting experimental challenge to observe the fast temperature decrease upon adiabatic demagnetization in these materials. Operation with relatively high magnetic fields makes rather questionable a technological application of the two gadolinium pyrochlores. The theoretical results suggest, though, a new direction for a search of suitable materials: a pyrochlore antiferromagnet with a factor of four smaller exchange constant than in Gd2Ti2O7 would be a more efficient refrigerant than Gd3Ga5O12, while operating in a similar range of fields and temperatures.

Real materials have, of course, more complicated magnetic interactions than just a nearest-neighbor exchange between classical spins. For example, in GGG a dipole-dipole interaction between neighboring spins amounts to 20–30% of the corresponding exchange energy.GGGphd  Long-range dipolar interactions play also a significant role in the two pyrochlore compounds.Raju ; Ramirez02  Influence of these and other extra perturbations on the magnetocaloric effect in frustrated magnets requires further theoretical analysis. Meantime, one can make several general remarks on the relative importance of various perturbations. The quantum effects omitted in the present study do not destroy the anomalous magnetothermodynamics near $`H_{\mathrm{sat}}`$. Similar to their classical counterparts quantum antiferromagnets on the three frustrated lattices also have flat branches of excitations. Honecker  The anomaly (7) is, nevertheless, affected by a weak dispersion ($`J_{\mathrm{disp}}`$) of soft excitations induced, e.g., by a second-neighbor exchange or dipole-dipole interactions. The magnetocaloric effect will saturate below $`T^{}J_{\mathrm{disp}}`$ at $`(S/H)_T1/\sqrt{T^{}}`$, which could be relatively large. Such a frustrated magnet can still be used, though less efficiently, as a refrigerant to reach temperatures below $`T^{}`$. A further analysis of a destructive role of various perturbations should focus, therefore, on an induced dispersion in flat magnon branches.

###### Acknowledgements.

I am greatly indebted to T. M. Rice, who’s questions have inspired this study. In the course of the work I benefited from fruitful discussions with A. I. Golov, C. L. Henley, A. Honecker, A. R. Muratov, and O. A. Petrenko.


format: text

text:
# On multisolitonic decay behavior of internal gravity waves

## 1 Introduction

There are papers (see S. B. Leble (1991) for a review) in which the problem of nonlinear waves evolution in stratified medium is in some way connected with soliton theory, see, e.g., Novikov et al (1984). We mean the conditions in which the physical (or environmental) conditions need account of nonlinearity and dispersion which joint action is of essential (sometimes - inevitable) importance. Evolution of internal gravity waves (IGW) in geophysical hydrodynamics and its further decay to turbulent spots is important example of such phenomena, Pedlosky (1979). The scenarios of this transition intensely discussed in literature of atmosphere and ocean physics Gavrilov, Fukao, Nakamura, Tsuda (1996). There are also publications about laboratory experiments devoted to IGW dynamics in tanks Segur, Finkel (1985). The problem have a link to a general problems of hydrodynamics as illustrated by Zeitunyan (1999). It is the Euler equations used as a basis of a description in which modes with different scales are generated within the process of a wave propagation over an initial density stratification taking into account its possible dynamic changes. This question relates to the behavior of a solution within the limiting procedure when a dispersion coefficient go to zero. In our case of internal waves it means the zero buoyancy frequency limit ($`N0`$ ) that in turn is proportional to the vertical density gradient. Let us mention also investigations of similar problems with a small dispersion parameter initiated by Lax and Levermore (1983) for the model of a single Korteweg-de Vries (KdV) equation. The dispersion in such model is described by the term with third derivative by a coordinate (x).

In this paper we study an evolution of a disturbance from an initial condition that model an excitation of the internal wave in the geometry and scales of McEwan (1983 a,b) experiments. We believe that the fine structures that appear in density distribution may be explained by a model based on a system of coupled Koreweg-de Vries (cKdV) equations. This system arises for coefficient functions related to a functional basis of transversal coordinate. It demonstrates the behaviors of solutions that have many features of multi-solitonic decay (see, again, Novikov et al (1984)) The wavetrains contain chains of such fluctuations of basic wave variables that propagate with velocity up to the most quick one. Our investigations show: interaction between the basic transverse guide modes does not destroy the general picture; only some phase changes appear. Of course the system of cKdV has higher order in time derivatives, the opposite directed wave train is launched, but the wave also has the multisoliton form.

The derivation of the system of cKdV may be considered as the direct application of Galerkin - Bubnov projection procedures that is successfully used in numerical simulation of initial-boundary problems with a good mathematical justification (the basis of this method are in Gottlieb and Orszag (1989), see, also Moser (1966)). After the variables division in the linearized problem some spectral problem appears. Its eigen functions of transverse variables form a basis. Hence the (infinite) 1+1 cKdV system is equivalent to the 2+1 basic Euler system under consideration within some class of initial conditions. While numeric simulation delivering we cut this infinite system restricting ourselves by some mode number, comparing few such examples. We, however believe that further extending of the system does not change the general behavior of the wavetrains. We plan checking of the hypothesis in the forthcoming paper but more concentrate now on the cKdV model that is also interesting by itself (e.g. Hirota R., Satsuma J. (1981) , A. Perelomova (2000)). It was studied mathematically via many approaches as described recently by Ayse Karasu (1997). There are integrable examples (R. Dodd, A.Fordy (1982)): it helps to test the numerical scheme to be introduce in this article.

For such systems we construct a finite-difference scheme (Sec. 3) and use it for numerical integration. The theorems about convergence and stability we formulate and prove here confirm the statements we made. The computations description we place in the section 4 and the proofs in Appendix.

In this article we do not touch the mixing problems in its full extent as well as interaction with the stratification mode (background or stationary state) itself. Such problems solution need account of boundary conditions and simulations for a longer time intervals. Account of such ingredients in the model could shade the main phenomena we intend to demonstrate. Therefore we concentrate attention of a reader on the initial stage of the wave field development. Already within this initial stage period the horizontal multi-solitonic fine structure prevails: it seems very important. Our general approach and the numerical scheme power allow to consider longer time to account reflections phenomena. A comparison of the evaluated stream function fields with ones of the McEwan experiments show many similar features, that, seems to verify the model from the point of the real fluid dynamics.

## 2 Formulation of the problem

### 2.1 Laboratory experiment of McEwan

Experiment were made using a rectangular plate-glass tank with dimensions shown in figure (1) below. This tank is filled by a linearly stratified salt solution. A wave-forcing dismountable paddle swung about horizontal lateral axis ia used .

Tank filled with a continuous stratified ambient state water

Fig. (1) A. D. McEwan experimental configuration.

Bring the fluid to a state of incipient breaking and photographing the evolution of these breaking events induced by internal wave action. McEwan gives his observation  along the following stages:
$`(a)`$ Overturning.
$`(b)`$ Development of interleaving microstructure.
$`(c)`$ Static stability is restored but microstructure is preserved.
$`(d)`$ Gravitation to an equilibrium has changed the surrounding density profile between extremum isopycnals.

We simulate only tile stage (b) above to show the obvious plots for the mixing phenomena.

## 3 Theoretical description and mathematical problem.

Consider the basic system of Euler equations for internal water waves in two dimensions (xz) with a stable stratified ambient state and the buoyancy frequency $`N(z)`$.

$$\begin{array}{cc}u_x+w_z=0& \\ \rho _ou_t=\rho _o(\underset{¯}{v},\underset{¯}{})up_x& \\ \rho _ow_t=\rho _o(\underset{¯}{v},\underset{¯}{})wp_z\rho ^{^{}}g& \\ T_t^{^{}}+w\overline{T}_z=(\underset{¯}{v},\underset{¯}{})T^{^{}}& \end{array}$$
(1)
where u, v are velocity components, $`\rho _o`$ is the density, $`p`$ is the pressure, $`\rho ^{^{}}g`$ is the body force due to stratifications, $`\overline{T}_z`$ is the vertical background temperature gradient and $`T^{^{}}`$ is the temperature variables.

In the commonly accepted approximations (incompressibility, Boussinesq approximation, etc.) the solution of the system (1) is constructed as the following representation for the current function$`\psi (x,z,t)`$, (See, for details Appendix A)

$$\psi (z,x,t)=\underset{n=1}{\overset{L}{}}Z^n(z)\theta ^n(x,t),$$
(2)
where $`Z^n(z)`$ are solutions of the correspondent Sturm-Liouville problem (Equation (19) in Appendix A)

$$Z_{zz}+\frac{N^2}{c_n^2}Z=0,Z(0)=Z(h)=0,$$
(3)
and describe a vertical shape of the wave modes. The linear propagation velocities, $`c_n`$, play the role of eigenvalues. The series (2) may be considered in the context of the Galerkin - Bubnov method. (for a basis of this method see, for example, Gottlieb and Orszag (1989)). The coefficient functions $`\theta ^n(x,t)`$ are solutions of the coupled KdV system (equation (27) in appendix A )

$$\theta _t^n+c_n\theta _x^n+\sigma \underset{m,k}{g_{m,k}^n}\theta ^m\theta _x^k+\beta ^2d_n\theta _{xxx}^n=0,$$
(4)
$`\sigma ,\beta ^2`$are scale parameters. Each solution $`\theta ^n`$ (x, t) stands for n-th mode wave amplitude as a function of the space $`x(\mathrm{},\mathrm{})`$ and time t coordinates. Here the constants $`g_{mkn},d_{nm}`$ are nonlinear and dispersion constants. The coefficients $`g_{mkn},d_{nm}`$ have been derived by Leble (1991) to have the formulae (equation (28) in appendix A )

$$g_{m,k}^n=\frac{\sigma N^2c_n^2}{2}\underset{h}{\overset{}{}}\left[\left(\frac{1}{c_m^2}+\frac{2}{c_k^2}\right)Z^kZ_z^m\frac{1}{c_mc_k}Z^mZ_z^k\right]z^ndz,d_n=\frac{c_n^3}{2N^2}$$
(5)
We specify the problem by zero boundary conditions for $`\psi `$ with respect to the variable $`z`$ as follows from (3). Physically it corresponds to the zero value of the vertical velocity components at the upper (z = h) and lower boundaries (z = 0) levels. We select a localized along x axis initial condition described by a smooth enough function that model the paddle motion as in the experiments of McEwan (1983a). The function is also chosen antisymmetric along $`x`$ in relation to the paddle axis centered in the middle of the tank. The time intervals of simulations are taken such that the initial disturbance decays essentially but does not reach the boundaries.

In the simplest model to be presented here we consider the case of constant buoyancy frequency. In the concrete numerical calculations the choice of the constants is the following $`N=1.23s^1`$ for a stratified salted water in a rectangular tank 50 cm long filled to a depth 25 cm (Ref. McEwan (1983a)).

## 4 Numerical method

The numerical scheme is based on the solution of the nonlinear finite (difference) system

$$c\left(\left(\theta ^n\right)_i^{j+\frac{1}{2}}\left(\theta ^n\right)_i^j\right)/\frac{\tau }{2}+c_n\left(\left(\theta ^n\right)_{i+1}^j\left(\theta ^n\right)_{i1}^j\right)/2h+\underset{k,m}{}g_{mkn}\left(\theta ^m\right)_i^j\left(\left(\theta ^k\right)_{i+1}^j\left(\theta ^k\right)_{i1}^j\right)/2h$$
(6)
$$+\left(d_nc_nh^2/6\right)\left(\left(\theta ^n\right)_{i+2}^j2\left(\theta ^n\right)_{i+1}^j+2\left(\theta ^n\right)_{i1}^j\left(\theta ^n\right)_{i2}^j\right)/2h^3=0,$$

where $`n`$, $`m`$, $`k`$ are the modes numbers; $`i`$ and $`j`$ are discrete space and time variables respectively. The time step is denoted by $`\tau `$ while the spatial step size by $`h`$. The equation (6) is accompanied with a discrete equation for the intermediate layer as:

$$c\left(\left(\theta ^n\right)_i^{j+1}\left(\theta ^n\right)_i^j\right)/\tau +c_n\left(\left(\theta ^n\right)_{i+1}^{j+\frac{1}{2}}\left(\theta ^n\right)_{i1}^{j+\frac{1}{2}}\right)/2h+\underset{k,m}{}g_{mkn}\left(\theta ^m\right)_i^{j+\frac{1}{2}}\left(\left(\theta ^k\right)_{i+1}^{j+\frac{1}{2}}\left(\theta ^k\right)_{i1}^{j+\frac{1}{2}}\right)/2h$$
(7)
$$+e_n\left(\left(\theta ^n\right)_{i+2}^{j+\frac{1}{2}}2\left(\theta ^n\right)_{i+1}^{j+\frac{1}{2}}+2\left(\theta ^n\right)_{i1}^{j+\frac{1}{2}}\left(\theta ^n\right)_{i2}^{j+\frac{1}{2}}\right)/2h^3=0$$

To support the results of further simulations to be presented here we would recover some mathematical ingredients of our work. Let $`(u^n)_i^j`$ $`=u^n(x_i,t^j)`$ denotes the grid function obtained from the exact solution $`u^n(x,t)`$ of the Cauchy problem for the cKdV system calculated at the grid points $`x_i`$, $`t^j`$. Then the discrete solution $`(\theta ^n)_i^j`$ of the system (6), (7 ) defines the discrepancy vector $`V^j`$ (residual function) with the components $`(v^n)_i^j=(\theta ^n)_i^j(u^n)_i^j.`$ Introduce the $`L_2`$ norm for $`V^j`$:

$$V^{j+1}=\left(\underset{i}{}\underset{n}{}\left[\left(v^n\right)_i^j\right]^2h\right)^{\frac{1}{2}}$$
(8)
The convergence of the scheme (6), (7 ) at $`(\tau ,h0)`$ follows from the

###### Theorem 1

Let the solution $`u^n(x,t)`$ of the Cauchy problem for the cKdV system with $`t[0,\mathrm{});x(\mathrm{},\mathrm{})`$ exists with the third continuous derivative with respect to time $`t`$, and the fifth derivative with respect to $`x`$. Then the following inequality takes place.

$$V^{j+1}\frac{1\sqrt{1\frac{4M\left|g_{m,k}^n\right|_{max}}{h^{\frac{3}{2}}}MO\left(\tau +h^2\right)}}{\frac{2M\left|g_{m,k}^n\right|_{max}}{h^{\frac{3}{2}}}}MO\left(\tau ^2+h^2\right),$$
(9)
where $`M`$ is a constant independent of $`\tau `$, $`h`$ and $`\left|g_{m,k}^n\right|_{max}=\underset{1n,m,kL}{\mathrm{max}}(\left|g_{m,k}^n\right|)`$; $`L`$ is quantity of modes in (LABEL:coupledkdv).

Therefore, the convergence follows from the inequality (9). The proof of the theorem is placed in the Appendix. We give the proof for more simple case then the only one of the numeric scheme equations (6) is used in simulations. The proof for the scheme (6)-(7) with a half-step is analogous but significantly more cumbersome, and we decided to give not it here, taking in mind publishing in the forthcoming special mathematical paper . This scheme and proof may be of interest even in the scalar (single KdV) case in spite of existence of numerous versions of such results. We believe the ours is effective as it follows from the tests and experiments presented here.

###### Remark 2

We have to note that evidently the existence of the fifth derivative by $`x`$ and the third derivative by $`\tau `$ are unnecessarily rigid conditions. These conditions are used for simplicity of consideration only, and may be easily weaken. One can replace them by the condition of the third order derivative by $`x`$ and the first order derivative by $`x`$ existence, then only the convergence is more slow. Moreover, one can accept a concept of weak solutions to be applied in applications. Usually such a concept allow to weaken the requirement of derivative existence on one order of derivatives; however the author have not verified it for the problem under consideration because necessity of such a delicate analysis is not evident.

## 5 Computation and results

When the coefficients in the Sturm-Liouville equation (3) are constant it has very simple general solution

$`Z^n=B_nSin\left(\frac{n\pi z}{h}\right),n=1,2,3,\mathrm{},L`$

which tends to zero at boundaries. The eigenvalues

$$c_n=\frac{Nh}{n\pi },n=1,2,3,\mathrm{}L.$$
(10)
have the sense of linear IGW velocities. Normalization is determined by $`\underset{0}{\overset{h}{}}\left(Z^n\right)^2N^2𝑑z=1`$ and gives

$`B_n=\left(\frac{2}{N^2h}\right)^{.5}`$, Hence

$$Z^n(z)=\left(\frac{2}{N^2h}\right)^{.5}Sin\left(\frac{n\pi z}{h}\right)$$
(11)
As it was discussed in the Sec.2 one can select the initial perturbation for the current function (1) which have the general form at the point $`t=0`$

$`\psi (z,x,0)=\underset{n=1}{\overset{L}{}}Z^n(z)\theta ^n(x,0)=\phi (x,z)=\phi _1(x)\phi _2(z).`$

Taking smooth, exponentially localized initial condition for $`\phi _1(x)`$ while $`\phi _2(z)`$ that model the impact of some rotating paddle motion as

$$\phi _1(x)=\frac{a}{Cosh(\frac{x}{l})}and\phi _2(z)=\left(\frac{2}{N^2h}\right)^{.5}Sech(b(zz_0))Tanh(b(zz_0))$$
(12)
The choice of the form and the constants a, b, l qualitatively reflects the paddle movement (we restrict the movement by some isolated pulse) and the numerical value for the amplitude is estimated also from the description of the experiment of McEwan (1983a). Then

$$(Z^j,\psi )=\underset{n=1}{\overset{L}{}}(Z^j,Z^n)\theta ^n(x,0)=(Z^j,\phi _2(z))\phi _1(x).$$
(13)
$$(Z^j,Z^n)=\underset{0}{\overset{h}{}}N^2Z^jZ^n𝑑z=\{1,(j=n)\},\{0,(jn)\}$$
(14)
$$(Z^j,\phi _2(z))=\underset{0}{\overset{h}{}}N^2Z^j\phi _2(z)𝑑z.$$
(15)
For $`h=0.25mandN=1.23s^1`$, taking the first five nonzero modes with coefficients from (15) and the numbers (2,4,6,8,10) while the odd numbers give zero, together with (14) and putting in (13) to obtain $`\theta (x,0)`$. The estimation of the modes energy satisfactory keep that of the initial perturbation. Then using (1) we get the initial perturbation $`\psi (z,x,0)`$. We give cross section plot at the middle of the tank Fig.(2) to show the antisymmetry and the percentage of energy content in the initial perturbation model

Fig.2 $`\phi _2(z)`$ initially and after Fourier representation

For the wave profile after time t we solve numerically system (4) with coefficients from (5), (10) to obtain $`\theta ^n(x,t)`$ then using (1) for $`\psi (z,x,t)`$ evaluation. For the chosen constants the calculated coefficients of nonlinear, dispersion terms and modes propagation velocities are picked together in the table

| $`g_{mk}^2`$ | | | | | $`g_{mk}^4`$ | | | | |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| m$`\backslash `$k | 2 | 4 | 6 | 8 | 10 |  | m$`\backslash `$k | 2 | 4 | 6 | 8 | 10 |  |
| 2 | 0 | 72.3 | 0 | 0 | 0 |  | 2 | 28.9 | 0 | 57.8 | 0 | 0 |  |
| 4 | 28.9 | 0 | 202.4 | 0 | 0 |  | 4 | 0 | 0 | 0 | 144.5 | 0 |  |
| 6 | 0 | 130 | 0 | 390.1 | 0 |  | 6 | 0 | 0 | 0 | 0 | 260 |  |
| 8 | 0 | 0 | 289 | 0 | 635.7 |  | 8 | 0 | 57.8 | 0 | 0 | 0 |  |
| 10 | 0 | 0 | 0 | 505.7 | 0 |  | 10 | 0 | 0 | 144.5 | 0 | 0 |  |

| $`g_{mk}^6`$ | | | | | $`g_{mk}^8`$ | | | | |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| m$`\backslash `$k | 2 | 4 | 6 | 8 | 10 |  | m$`\backslash `$k | 2 | 4 | 6 | 8 | 10 |
| 2 | 0 | 33.7 | 0 | 53 | 0 |  | 2 | 0 | 0 | 36.1 | 0 | 50.6 |
| 4 | 48.2 | 0 | 0 | 0 | 125.2 |  | 4 | 0 | 57.8 | 0 | 0 | 0 |
| 6 | 0 | 0 | 0 | 0 | 0 |  | 6 | 65 | 0 | 0 | 0 | 0 |
| 8 | -19.3 | 0 | 0 | 0 | 0 |  | 8 | 0 | 0 | 0 | 0 | 0 |
| 10 | 0 | 24 | 0 | 0 | 0 |  | 10 | -36 | 0 | 0 | 0 | 0 |

| $`g_{mk}^{10}`$ | | | | |  |
| --- | --- | --- | --- | --- | --- |
| m$`\backslash `$k | 2 | 4 | 6 | 8 | 10 |
| 2 | 0 | 0 | 0 | 37.6 | 0 |
| 4 | 0 | 0 | 63.6 | 0 | 0 |
| 6 | 0 | 78 | 0 | 0 | 0 |
| 8 | 80.9 | 0 | 0 | 0 | 0 |
| 10 | 0 | 0 | 0 | 0 | 0 |

$`d_{22}`$ =$`.00004`$, $`d_{44}`$ =$`.000005`$, $`d_{66}`$ =$`.000002`$, $`d_{88}`$=$`.000001`$, $`d_{1010}`$=$`.0000003`$,

$`c_2`$=$`.05`$, $`c_4`$=$`.025`$ , $`c_6`$=$`.016`$, $`c_8`$=$`.012`$, $`c_{10}`$=$`.0098`$

For the case of linearized equations the wave profile after time t=0.02 s may be obtained directly from the same system if one put nonlinear constants to zero. The evaluation result is shown in the following contour plot Fig.(3) for the upper half of the tank.

Fig.3 The wave profile 3-dimensional plot in linear case with dispersion for the upper half of the tank (12.5 cm in z and 50 cm in x directions)

For the nonlinear case we show one-dimensional (x) profiles Fig.(4.a,b,c,d,e) for the second, forth, sixth, eighth and tenth modes respectively

Fig.(4.a) Second mode, initially (dotted) and at t=0.02 sec(continuous).

Fig.(4.b) Fourth mode, initially (dotted) and at t=0.02 sec(continuous).

Fig.(4.c) Sixth mode, initially (dotted) and at t=0.02 sec(continuous).

Fig.(4.d) Eighth mode, initially (dotted) and at t=0.02 sec(continuous).

Fig.(4.e) Tenth mode, initially (dotted) and at t=0.02 sec(continuous).

Fig.4 One-dimensional (x) plots of the wave profile in nonlinear case.

The contour plot for the wave profile at $`t=0.02`$ sec is presented below (Fig.5) for the same half of the tank.

Fig.4 Three dimensional (x,y,z) plot of the wave profile for the upper half of the tank.

The horizontal numbering indicates the number of points used in plot while the vertical

show the actual z profile of the waves

All the plots for the modes and for the sum (current function) show the behavior that is typical for the multisolitonic perturbation. It looks likely a decay of the initial condition to solitons in the single KdV equation theory. The difference between linear dispersion and nonlinear looks crucial. Even within the short time interval we present here the horizontal scale of the disturbance radically changes. The final scale is defined by constants of the cKdV system and the amplitude of initial condition. The process of the wave propagation is accompanied by interaction that imply the energy transfer between modes. This phenomenon may be considered as a possible reason for the vertical fine structure generation (Ref. Leble (1991) ) The combination of the fine structures may explain the experiments of McEwan (1983 a,b).

## 6 Conclusion

We conclude with the following remarks about restrictive assumptions we made and perspectives. First in the framework of clarity of the plotting and interpretation we considered only one direction of propagating waves. We would repeat that one can elongate the time of simulations to have more interactions but such picture looks cumbersome.

The relatively small number of modes is taken into account because it conserves most the energy. The number may be increased but at the times we consider the effect of small-vertical-scales modes generation may be neglected. Returning to the general significance in the theory of the Euler equations validity mentioned in the introduction, one can easily put the important question. The question arises when we analyze the expressions for the constants in the of equations cKdV system. In a limit of the homogeneous medium, or $`N0`$, the nonlinear constant increase while the dispersion and velocities go to zero. It means the growth of solitons number and more strong interaction between modes. The discussion of this phenomenon and ideas of adequate regularization will be considered in the forthcoming publication.

## 7 Appendix A <br>Theoretical description and mathematical problem.

Consider Euler equations that describe IGW in two dimensions (x,z) stratified water

$$\begin{array}{cc}u_x+w_z=0& \\ \rho _ou_t=\rho _o(\underset{¯}{v},\underset{¯}{})up_x& \\ \rho _ow_t=\rho _o(\underset{¯}{v},\underset{¯}{})wp_z\rho ^{^{}}g& \\ T_t^{^{}}+w\overline{T}_z=(\underset{¯}{v},\underset{¯}{})T^{^{}}& \end{array}$$
(16)
where u, v are velocity components, $`\rho _o`$ is the density, $`p`$ is the pressure, $`\rho ^{^{}}g`$ is the body force due to stratifications, $`\overline{T}_z`$ is the vertical background temperature gradient and $`T^{^{}}`$ is the temperature variables.

Combining equations in (16) and using the state equation for liquid $`\rho ^{^{}}=\rho _o\alpha T^{^{}},`$ $`\alpha `$ $`=\frac{\rho _z}{\rho \overline{T}_z}`$is the coefficients of thermal expansion $`\overline{T}_z=\frac{N^2}{\alpha g}`$we obtain

$$\mathrm{\Delta }w_{tt}+N^2w_{xx}N^2\left[(\underset{¯}{v},\underset{¯}{})_0^tw𝑑t\right]_{xx}+\left[(\underset{¯}{v},\underset{¯}{})w_z𝑑x\right]_{txz}+\left[(\underset{¯}{v},\underset{¯}{})w\right]_{txx}=0$$
(17)
Going to the dimensionless variables by the following formulas

$`x_i=\lambda _ix_i^{},t=\frac{2\pi }{\overline{N}\beta }t^{},u=\frac{\lambda _z\overline{N}}{2\pi }u^{},w=\frac{\beta \lambda _z\overline{N}}{2\pi }w^{},\alpha T^{}=T,N=\frac{\overline{N}}{2\pi }N^{}`$

So $`\left(w/x_i\right)=\left(1/\lambda _i\right)\left(w/x_i^{}\right),\left(w/t\right)=\left(\overline{N}\beta /2\pi \right)\left(w/t^{^{}}\right)`$

Substitute in (17) and omit primes for simplicity

$`\lambda _z\left(\frac{\overline{N}\beta }{2\pi }\right)^3\left(\frac{w_{xx}}{\lambda _x^2}+\frac{w_{zz}}{\lambda _z^2}\right)_{tt}+\left(\frac{\overline{N}}{2\pi }\right)^3\left(\beta N\right)^2\left(\frac{w_{xx}}{\lambda _x}\right)\left(\frac{\overline{N}}{2\pi }\right)^3\frac{\left(\beta N\right)^2}{\lambda _x}\left[u_0^tw_x𝑑t+w_0^tw_z𝑑t\right]_{xx}`$

$`+\left(\frac{\overline{N}\beta }{2\pi }\right)^3\left\{\frac{\lambda _z}{\lambda _x^2}\left[uw_x+ww_z\right]_{xxt}+\frac{1}{\lambda _z}\left[uw_z+w\left(w_z𝑑x\right)_z\right]_{xzt}\right\}=0`$

Introducing the current function , $`w=\sigma \psi _x,u=\sigma \psi _z`$ , Integrate w.r.t. x

$$\psi _{zztt}+N^2\psi _{xx}=\beta ^2\psi _{xxtt}\sigma \left[\psi _z\psi _{xz}\psi _x\psi _{zz}\right]_{zt}+\sigma N^2\left[\psi _z_0^t\psi _{xx}𝑑t\psi _x_0^t\psi _{xz}𝑑t\right]_x$$
(18)
Substitute in (18) by $`\psi =\underset{𝑚}{}z^m\theta ^m,`$multiply by$`z^n`$ , inegrate w.r.t. z and using The separation of variables that gives

$$z_{zz}^n=\frac{N^2}{c_n^2}z^n$$
(19)
So we have

$$\theta _{tt}^nc_n^2\theta _{xx}^n=\frac{c_n^2\beta ^2}{N^2}\theta _{xxtt}^n+\sigma c_n^2\underset{m,k}{}a_{m,k}^n\left(\theta ^m\theta _x^k\right)_t+\left[b_{m,k}^n\theta ^m\theta _t^k+e_{m,k}^n\theta _x^m\underset{0}{\overset{𝑡}{}}\theta _x^kdt\right]_x$$
(20)
where$`a_{m,k}^n=N^2\underset{h}{\overset{}{}}\left(\frac{1}{c_k^2}+\frac{1}{c_m^2}\right)z_z^mz^kz^ndz,b_{m,k}^n=\frac{N^2}{c_k^2}\underset{h}{\overset{}{}}z_z^mz^kz^ndz,e_{m,k}^n=N^2\underset{h}{\overset{}{}}z^mz_z^kz^ndz`$

The equations of the separated propagated modes are obtained as follow
Let $`\theta _t^n=u^n,c_n\theta _x^n=v^n`$

So (20) becomes

$$\begin{array}{cc}u_t^nc_nu_x^n=\frac{c_n^3\beta ^2}{N^2}v_{xxx}+Nonlinearterms& \\ v_t^nc_nu_x^n=0& \end{array}$$
(21)
Using the projection operator P, $`P_+=\frac{1}{2}\left(\begin{array}{c}\mathrm{1\hspace{0.17em}\; \hspace{0.17em}\; 1}+\frac{c^2\beta ^2}{2N^2}_x^2\\ 1\frac{c^2\beta ^2}{2N^2}_x^2\mathrm{\hspace{0.17em}\; 1}\end{array}\right),P_{}=\frac{1}{2}\left(\begin{array}{c}11\frac{c^2\beta ^2}{2N^2}_x^2\\ 1+\frac{c^2\beta ^2}{2N^2}_x^2\mathrm{\hspace{0.17em}\; 1}\end{array}\right)`$

$$\begin{array}{cc}P_+\left(\begin{array}{c}u^n\\ v\end{array}\right)=\left(\begin{array}{c}\phi ^{n+}\\ k\phi ^{n+}\end{array}\right),P_{}\left(\begin{array}{c}u\\ v\end{array}\right)=\left(\begin{array}{c}\phi ^{}\\ k\phi ^{}\end{array}\right)& \\ u=\phi ^{n+}+\phi ^n,v=\left(\phi ^{n+}\phi ^n\right)\frac{c^2\beta ^2}{2N^2}_x^2\left(\phi ^{n+}\phi ^n\right)& \end{array}$$
(22)
Operating $`P_+,P_{}`$on (21) and using (22) we obtain the equations for the separated modes $`\phi ^{n+},\phi ^n`$

$$\begin{array}{cc}\phi _t^{n+}c_n\phi _x^{n+}\frac{c^3\beta ^2}{2N^2}\phi _{xxx}^{n+}\frac{\sigma c_n^2}{2}\underset{m,k}{}a_{m,k}^n((\phi ^++\phi ^{})^mdt.\frac{\left(\phi ^+\phi ^{}\right)^k}{c_k})_t& \\ +[b_{m,k}^n(\phi ^++\phi ^{})^mdt.(\phi ^++\phi ^{})^k+e_{m,k}^n\frac{\left(\phi ^+\phi ^{}\right)^m}{c_m}\underset{0}{\overset{𝑡}{}}\frac{\left(\phi ^+\phi ^{}\right)^k}{c_k}dt]_x& \end{array}$$
(23)
$$\begin{array}{cc}\phi _t^n+c_n\phi _x^n+\frac{c^3\beta ^2}{2N^2}\phi _{xxx}^n\frac{\sigma c_n^2}{2}\underset{m,k}{}a_{m,k}^n((\phi ^{n+}+\phi ^n)^mdt.\frac{\left(\phi ^+\phi ^{}\right)^k}{c_k})_t& \\ +[b_{m,k}^n(\phi ^++\phi ^{})^mdt.(\phi ^++\phi ^{})^k+e_{m,k}^n\frac{\left(\phi ^{n+}\phi ^n\right)^m}{c_m}\underset{0}{\overset{𝑡}{}}\frac{\left(\phi ^+\phi ^{}\right)^k}{c_k}dt]_x& \end{array}$$
(24)
Let $`\phi ^{n+}=\mathrm{\Psi }_t^{n+},\phi ^n=\mathrm{\Psi }_t^n`$ So (23), (24) becomes

$$\begin{array}{cc}\mathrm{\Psi }_t^{n+}c_n\mathrm{\Psi }_x^{n+}=\frac{c^3\beta ^2}{2N^2}\mathrm{\Psi }_{xxx}^{n+}+\frac{\sigma c_n^2}{2}\underset{m,k}{\frac{a_{m,k}^n}{c_k}}\left(\mathrm{\Psi }^++\mathrm{\Psi }^{}\right)^m\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)_x^k& \\ +b_{m,k}^n\left(\mathrm{\Psi }^++\mathrm{\Psi }^{}\right)^m.\left(\phi ^++\phi ^{}\right)_x^k+\frac{e_{m,k}^n}{c_kc_m}\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)_x^m\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)^k& \end{array}$$
(25)
$$\begin{array}{cc}\mathrm{\Psi }_t^{n+}+c_n\mathrm{\Psi }_x^{n+}=\frac{c^3\beta ^2}{2N^2}\mathrm{\Psi }_{xxx}^{n+}+\frac{\sigma c_n^2}{2}\underset{m,k}{\frac{a_{m,k}^n}{c_k}}\left(\mathrm{\Psi }^++\mathrm{\Psi }^{}\right)^m\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)_x^k& \\ +b_{m,k}^n\left(\mathrm{\Psi }^++\mathrm{\Psi }^{}\right)^m.\left(\phi ^++\phi ^{}\right)_x^k+\frac{e_{m,k}^n}{c_kc_m}\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)_x^m\left(\mathrm{\Psi }^+\mathrm{\Psi }^{}\right)^k& \end{array}$$
(26)
In next work we will consider both directed modes. In this work we evaluate only for short time to describe the phenomena. so we consider only one direction that have the form

$$\mathrm{\Psi }_t^n+c_n\mathrm{\Psi }_x^n+\sigma \underset{m,k}{g_{m,k}^n}\mathrm{\Psi }^m\mathrm{\Psi }_x^k+\beta ^2d_n\mathrm{\Psi }_{xxx}^n=0,$$
(27)
$`\sigma ,\beta ^2`$are scale parameters ,

$$a_{m,k}^n=\frac{\sigma N^2c_n^2}{2}\underset{h}{\overset{}{}}\left[\left(\frac{1}{c_m^2}+\frac{2}{c_k^2}\right)Z^kZ_z^m\frac{1}{c_mc_k}Z^mZ_z^k\right]z^ndz,d_n=\frac{c_n^3}{2N^2}$$
(28)
## 8 Appendix B <br>Proof of the convergence.

Now we realize that the finite-difference method used gives us a correct solution of the problem if the time and space steps is enough small and we show how we can easily estimate an error of the numerical method. At first step we prove that a solution $`(\theta _i^n)^j`$ of finite-difference equations (6), 7) converges to some solution of equation (LABEL:coupledkdv) as $`\tau ,h0`$.

We give here a proof for more simple case for brevity when the scheme is based on equation (6). The proof for equations (6), 7) is analogous but more cumbersome. Therefore, now we consider the difference equations

$$\frac{\left(\theta ^n\right)_i^{j+1}\left(\theta ^n\right)_i^j}{\tau }+c_n\frac{\left(\theta ^n\right)_{i+1}^j\left(\theta ^n\right)_{i1}^j}{2h}+\underset{m,k}{}g_{m,k}^n\left(\theta ^m\right)_i^j\frac{\left(\theta ^k\right)_{i+1}^j\left(\theta ^k\right)_{i1}^j}{2h}$$
(29)
$$+e_n\frac{\left(\theta ^n\right)_{i+2}^j2\left(\theta ^n\right)_{i+1}^j+2\left(\theta ^n\right)_{i1}^j\left(\theta ^n\right)_{i2}^j}{2h^3}=0$$

Here $`h`$ is a space grid step, $`i`$ is a discrete space variable, $`j`$ is a discrete time, $`\tau `$ is a time step, $`n`$ is a mode number. We suppose that an exact solution is a smooth one and therefore

$`{\displaystyle \frac{(\theta _i^n)^{j+1}(\theta _i^n)^j}{\tau }}=\theta _t^n|_{x=ih}+O(\tau ){\displaystyle \frac{(\theta _{i+1}^n)^j(\theta _{i1}^n)^j}{2h}}=\theta _x^n|_{x=ih}+O(h^2)`$ (30)
$`{\displaystyle \frac{(\theta _{i+2}^n)^j2(\theta _{i+1}^n)^j+2(\theta _{i1}^n)^j(\theta _{i2}^n)^j}{2h^3}}=\theta _{xxx}^n|_{x=ih}+O(h^2)`$
In this way, equations (29) approximates equations (LABEL:coupledkdv) with an error of the order of $`O(\tau +h^2)`$.

We know that a single KdV equation $`u_t+uu_x+u_{xxx}=0`$ has the conservation law $`_{\mathrm{}}^{\mathrm{}}u^2(x,t)𝑑x=const`$. So we can consider a solution of a single KdV equation as an element of $`L_2`$-space, parameterized by time $`t`$. Equations (LABEL:coupledkdv) are of the same structure as a simple KdV equation, so it is naturally to develop the proof in the functional space with $`L_2`$-norm.

Let $`\left(\theta ^n\right)^j`$ denote a column

$$\left(\theta ^n\right)^j\left(\begin{array}{c}\mathrm{}\\ \left(\theta ^n\right)_{i2}^j\\ \left(\theta ^n\right)_{i1}^j\\ \left(\theta ^n\right)_i^j\\ \left(\theta ^n\right)_{i+1}^j\\ \left(\theta ^n\right)_{i+2}^j\\ \mathrm{}\end{array}\right),n=1,2,\mathrm{},l$$

and let $`\mathrm{\Theta }^j`$ denote a column consisting of columns $`\left(\theta ^n\right)^j`$ introduced above:

$$\mathrm{\Theta }^j\left(\begin{array}{c}\left(\theta ^1\right)^j\\ \left(\theta ^2\right)^j\\ \left(\theta ^3\right)^j\\ \mathrm{}\\ \left(\theta ^l\right)^j\end{array}\right)$$

Let $`u^n(x,t)`$ be an exact solution of (LABEL:coupledkdv). Then $`(u_i^n)^ju^n(x_i=ih,t=j\tau )`$ is an exact solution in the point $`x_i`$ of the grid at the moment $`j\tau `$. Let

$$\left(u^n\right)^j\left(\begin{array}{c}\mathrm{}\\ \left(u^n\right)_{i2}^j\\ \left(u^n\right)_{i1}^j\\ \left(u^n\right)_i^j\\ \left(u^n\right)_{i+1}^j\\ \left(u^n\right)_{i+2}^j\\ \mathrm{}\end{array}\right),n=1,2,\mathrm{},l,\text{and}U^j\left(\begin{array}{c}\left(u^1\right)^j\\ \left(u^2\right)^j\\ \left(u^3\right)^j\\ \mathrm{}\\ \left(u^l\right)^j\end{array}\right)$$

The difference $`V^j=\mathrm{\Theta }^jU^j`$ is an error of the numerical method. The $`L_2`$-norm of $`V^j`$ is

$$V^j=\left(\underset{n=1}{\overset{l}{}}\underset{i=\mathrm{}}{\overset{+\mathrm{}}{}}\left[\left((v^n)_i^j\right)^2h\right]\right)^{\frac{1}{2}}.$$
(31)
The difference solution $`\left(\theta ^n\right)_i^j`$ converges to the exact solution $`\left(u^n\right)_i^j`$ if $`V^j0`$ as $`\tau ,h0`$.

We substitute (LABEL:App-Error) into (29) and obtain an equation for the error $`(v^n)_i^j`$

$$\frac{\left(v^n\right)_i^{j+1}\left(v^n\right)_i^j}{\tau }+c_n\frac{\left(v^n\right)_{i+1}^j\left(v^n\right)_{i1}^j}{2h}+\underset{m,k}{}(g_{m,k}^n\left(u^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h}+g_{m,k}^n\left(v^m\right)_i^j\frac{\left(u^k\right)_{i+1}^j\left(u^k\right)_{i1}^j}{2h}$$
(32)
$$+g_{m,k}^n\left(v^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h})+e_n\frac{\left(v^n\right)_{i+2}^j2\left(v^n\right)_{i+1}^j+2\left(v^n\right)_{i1}^j\left(v^n\right)_{i2}^j}{2h^3}=(\frac{\left(u^n\right)_i^{j+1}\left(u^n\right)_i^j}{\tau }+$$

$$c_n\frac{\left(u^n\right)_{i+1}^j\left(u^n\right)_{i1}^j}{2h}+\underset{m,k}{}g_{m,k}^n\left(u^m\right)_i^j\frac{\left(u^k\right)_{i+1}^j\left(u^k\right)_{i1}^j}{2h}+e_n\frac{\left(u^n\right)_{i+2}^j2\left(u^n\right)_{i+1}^j+2\left(u^n\right)_{i1}^j\left(u^n\right)_{i2}^j}{2h^3})$$

Pick out a linear part of expression (32) and introduce for convenience an operator $`T^j`$ by the expression

$$\left(v^n\right)_i^j\tau (c_n\frac{\left(v^n\right)_{i+1}^j\left(v^n\right)_{i1}^j}{2h}\underset{m,k}{}(g_{m,k}^n\left(u^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h}+g_{m,k}^n\left(v^m\right)_i^j\frac{\left(u^k\right)_{i+1}^j\left(u^k\right)_{i1}^j}{2h})+$$
(33)
$$e_n\frac{\left(v^n\right)_{i+2}^j2\left(v^n\right)_{i+1}^j+2\left(v^n\right)_{i1}^j\left(v^n\right)_{i2}^j}{2h^3})=_{r=\mathrm{}}^{\mathrm{}}(T^{j+1}\left)_{ir}\right(v^n)_r^j,n=1,2,3,\mathrm{}l.$$

Because $`\left(u^n\right)_i^j`$ is an exact solution of differential equations, the right part of (32) is a small of the order $`O\left(\tau +h^2\right).`$ Therefore, we can rewrite (32) in the following way

$$\left(v^n\right)_i^{j+1}\underset{r}{}\left(T^{j+1}\right)_{ir}\left(v^n\right)_i^j+\tau \underset{m,k}{}g_{m,k}^n\left(v^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h}=\tau \left(f^n\right)_i^j,$$

$$\left(f^n\right)_i^j=\underset{m,k=1}{\overset{l}{}}g_{m,k}^n\left(v^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h}+O\left(\tau +h^2\right)$$

We need an estimate of a norm of $`T^j`$. It is valid

###### Lemma 3

If $`u^n(x,t)`$ is a function with a bounded derivative $`\left|\frac{}{x}u^n(x,t)\right|`$, then for any $`b>0`$ and for all

$$\tau bh^6,$$
(34)
there exists a constant $`A`$, independent of $`\tau `$, $`h`$, $`j`$, such that $`T^j_s\mathrm{exp}(A\tau )`$.

Here $`T^j_s`$ is a spectral norm of the operator $`T^j`$, i.e., the norm of $`T^j`$ that is induced by $`L_2`$-norm $`V^j`$. We will omit the subscribe ”$`s`$” in the following and will simply write $`T^j`$. The proof of Lemma 3 is a standard one, but cumbersome. Therefore we do not give the proof of Lemma 3 here and recommend reader books , .

Let $`F^j`$ be the column constructed from $`(f_i^n)^j`$ in the same manner as $`\mathrm{\Theta }^j`$ constructed from $`(\theta _i^n)^j`$. We need a norm of $`F^j`$; simple estimation gives

$$F^j=\left(\underset{n=1}{\overset{l}{}}\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}\left[\left(f^n\right)_i^j\right]^2h\right)^{\frac{1}{2}}=G\left(\underset{n=1}{\overset{l}{}}\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}\left[\underset{m,k=1}{\overset{l}{}}\left(v^m\right)_i^j\frac{\left(v^k\right)_{i+1}^j\left(v^k\right)_{i1}^j}{2h}\right]^2h\right)^{\frac{1}{2}}+O(\tau +h^2)$$
(35)
$$G\underset{n=1}{\overset{l}{}}\left(\underset{m=1}{\overset{l}{}}\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}\left[\left(v^m\right)_i^j\right]^2h\underset{k=1}{\overset{l}{}}\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}\left[\left(v^k\right)_i^j\right]^2h\right)^{\frac{1}{2}}\frac{1}{h^{\frac{3}{2}}}+O\left(\tau +h^2\right)\frac{Gl}{h^{\frac{3}{2}}}\left(V^j\right)^2+O\left(\tau +h^2\right)$$

where $`G=\underset{1n,m,kl}{\mathrm{max}}\left|g_{m,k}^n\right|`$.

So, using the operator $`T^j`$, we can rewrite (32) as

$$V^{j+1}=T^{j+1}V^j+\tau F^j.$$

The following estimate is valid

$$V^{j+1}T^{j+1}V^j+\tau F^j$$
(36)
If we consequently substitute (36) into itself, we get

$`V^{j+1}T^{j+1}V^j+\tau F^jT^{j+1}T^jV^{j1}+\tau (T^{j+1}F^{j1}+F^j)`$ (37)
$`T^{j+1}T^jT^{j1}V^{j2}+\tau (T^{j+1}T^jF^{j2}+T^{j+1}F^{j1}+F^j`$
$`e^{A\tau j}V^0+\tau (e^{A\tau (j1)}F^0+e^{A\tau (j2)}F^1+\mathrm{}F^j)`$
$`e^{A\tau j}V^0+M\underset{kj}{\mathrm{max}}(F^k)`$
$`e^{A\tau j}V^0+MF^{j+1},`$
$`M=\tau {\displaystyle \frac{e^{At_{\mathrm{max}}}1}{e^{A\tau }1}}.`$
We have used a known formula for the sum of a geometric series. Here $`t_{\mathrm{max}}`$ is the time of simulation: $`0tt_{\mathrm{max}}`$ .

Substitution of (35) instead of $`F^j`$ into (37) yields the inequality

$$V^{j+1}e^{A\tau j}V^0+M\left(\frac{Gl}{h^{\frac{3}{2}}}\left(V^{j+1}\right)^2\right)+O\left(\tau +h^2\right)$$

Because of the initial conditions we know exactly, we have $`(v_i^n)^0=0`$. So, taking into consideration $`V^0=0`$ we obtain the solution of the inequality

$$V^{j+1}\frac{1\sqrt{1\frac{4MGl}{h^{\frac{3}{2}}}MO\left(\tau +h^2\right)}}{\frac{2MGl}{h^{\frac{3}{2}}}}MO\left(\tau +h^2\right)$$

i.e. $`V^{j+1}0`$ as $`\tau `$, $`h0`$. The convergence is proved.

The error of the method considered is of order $`O\left(\tau +h^2\right)`$, if the solution under consideration is enough smooth one. The method converges for any $`\tau ,h0`$, but the rate of converging is dependent of relation between $`\tau `$ and $`h`$: the more $`\frac{\tau }{h^6}`$ the better. If the solution is not smooth, then we must consider the solution in the sense of a distribution theory; this case is not study yet.

The proof of the convergence for the scheme (6), (7 ) is analogous, but more cumbersome and tedious. The error of the method (6), (7 ) is $`O\left(\tau ^2+h^2\right)`$, and instead of condition (34 ) the condition $`\tau bh^4`$ is used. The main idea of the proof is the same one: at small $`\tau `$, $`h`$ the error due to difference approximation of the nonlinear terms influence on the result more less then the errors due to difference approximation of dispersion terms; and all difficult estimates only are necessary to show this fact. Therefore, choice of the optimal step $`tau`$ may be accomplished with the help of consideration of simple difference equation with dispersion term only. We had used this fact when made numerical experiments.


text:
# 1 Introduction

## 1 Introduction

It is generally believed that the precision electroweak(EW) data are in good agreement with the standard model(SM) of particle physics. Yet one has to admit that during the last couple of years a subset of the data has made the situation somewhat uncomfortable for the SM. If only a small subset of a large volume of data is in strong disagreement with theory, such that its removal improves the quality of the fit, then one would normally believe that either statistical fluctuations or some hitherto unknown systematic errors might have affected this subset. A conservative approach would then be to keep the subset in the cold storage awaiting improved statistics or a better understanding of the systematics, rather than invoking a new physics model.

The situation, unfortunately, is not that simple in the present case. It is now well known that the effective leptonic weak mixing angle ( given by $`sin^2\theta _W^l`$ = $`x_W^l`$) determined from the leptonic asymmetries measured by the LEP groups and SLD are in severe disagreement with the same parameter as determined by the hadronic asymmetries. The discrepancy in the global analysis becomes more prominent if the effective on-shell weak mixing angle extracted from recent neutrino scattering data is taken into account. In Ref. several standard model fits were performed (see Table 13.2 ) and it was noted that the $`\chi ^2`$/d.o.f was large in each case. The large dispersion in the fitted value of $`x_W^l`$ from various asymmetries, responsible for the poor fit, was interpreted as the result of fluctuations in one or more of the experimental inputs and the issue was not pursued any further. Unfortunately no improvement in statistics is expected in the near future. Most of the analyses based on the measurements at the $`Z`$-pole are now more or less final. With more than 700 pb-1 of data per experiment at energies higher than the $`Z`$-pole an improved accuracy in $`m_W`$ measurement is expected which in turn may lead to a better understanding of the SM vis a vis LEP data.

If the hadronic asymmetries are discarded the quality of the fit improves dramatically, as expected, but the mass of the Higgs boson obtained from the fit turns out to be much smaller than the lower limit obtained from the direct searches at LEP. Given the theoretical uncertainties (unknown higher order corrections, the precise value of an important input - the fine structure constant evaluated at $`m_Z`$ etc.) nicely reviewed in Ref., the possibility of statistical fluctuations and poorly understood systematics, it is not impossible that the data can still be in agreement with the SM. However, as analyzed in great details in Ref., something totally unexpected has to happen. For example, if statistical fluctuation is the possible explanation then it is imperative that not only the measurements in disagreement ( some of the hadronic asymmetries, say ) with the SM but also the ones which have been thought to be the evidence for the SM for so many years must involve large fluctuations. It was, therefore, argued in Ref. that new physics seems to be the favoured solution, although the evidence is not fully conclusive.

Several new physics models have already been proposed as possible solutions of the alleged conflict between the SM and the data (see Ref. for a list of further references). In this paper the focus will be on the supersymmetric solution. This solution seems to be attractive because one does not need supersymmetry only to ameliorate the malady in the precision electroweak data. It is needed to answer deeper issues. It solves the hierarchy puzzle that haunts non-supersymmetric grand unified theories and facilitates the coupling constant unification. The supersymmetric extension of the minimal standard model (MSSM), therefore, seems to be a well motivated step beyond the SM. It was shown in Ref. that the quality of the fit to the precision electroweak data as well as the agreement with the lower bound on the Higgs mass improve in the MSSM with light slepton-sneutrinos. In particular, the MSSM with sneutrino ($`\stackrel{~}{\nu }`$)s having mass in the range $`5580`$ GeV seems to be preferred by the data.The left slepton masses ($`m_{\stackrel{~}{\mathrm{}}_L}`$) are related to $`m_{\stackrel{~}{\nu }}`$ by the SU(2) breaking D-term in a model independent way. In order to make them heavier than the LEP bound moderate or large $`\mathrm{tan}\beta `$ was needed The fit favoured much heavier squarks (mass $``$ 1 TeV), was practiclly independent of $`m_A`$ (the mass of the pseudoscalar Higgs boson) and $`m_{\stackrel{~}{\mathrm{}}_R}`$ (the mass of the right slepton). Gaugino mass unification was assumed although the U(1) gaugino mass parameter had little impact. Mass of the lighter chargino ($`m_{\stackrel{~}{\chi }^\pm })<`$ 150 GeV improved the fit although higher masses also had reasonable agreement with the data. Left sleptons and sneutrinos belonging to different generations were assumed to be degenerate. .

The emphasis of this paper would be on the possibility of testing the light sneutrino hypothesis via direct searches using existing data. This is especially important since any immediate improvement in the indirect test of this scenario using electroweak data is not feasible. It is gratifying to note that the existing LEP data on single $`\gamma `$ \+ missing energy events  can indeed shed light on this issue.

However, before taking up the main issue, we want to review briefly a related topic. Is there a theoretically well motivated supersymmetry breaking mechanism which leads naturally to the light slepton scenario? It was noted in Ref. that this scenario cannot be accommodated in the popular minimal supergravity (mSUGRA) model where supersymmetry breaking is driven by gravity mediated interactions leading to a common scalar mass ($`m_0`$) at a scale which is often assumed, somewhat arbitrarily, to be the GUT scale ($`M_G`$). In such models the masses of the sleptons and the sneutrinos are correlated. As a result $`m_{\stackrel{~}{\nu }}`$ in the above mass range would inevitably lead to a somewhat lighter right slepton with mass ($`m_{\stackrel{~}{\mathrm{}}_R}`$) in conflict with the existing LEP lower bound.

This observation has recently been quantified by the ALEPH collaboration. Their work is based on a model, similar to mSUGRA, with a common $`m_0`$ and a common gaugino mass at the GUT scale. The limit is based on the data for various sfermion-gaugino searches. No direct sneutrino signal was searched for. Yet from the absence of any slepton signal they obtained a model dependent lower bound on the sneutrino mass ($`m_{\stackrel{~}{\nu }}>`$ 84 GeV) by exploiting the correlation between slepton and sneutrino masses.

Possible alternative scenarios with light sleptons were qualitatively discussed in Ref.. It was noted that if $`m_0`$ is generated at the Planck scale ($`M_P`$) instead of $`M_G`$ then the running between $`M_P`$ and $`M_G`$ may indeed lead to a somewhat larger $`m_{\stackrel{~}{\mathrm{}}_R}`$ at the weak scale. This avoids the conflict with the LEP bound. In fact an inspection of this running within the framework of an SU(5) SUSY GUT, as given in Ref., would encourage this scenario. The other viable model proposed was the anomaly mediated supersymmetry breaking (AMSB) model. In the latter model, without any additional assumption, the slepton-sneutrino masses turn out to be tachyonic. In the simplest version of this model one adds a common soft breaking term making the mass squared terms positive. The slepton mass can, therefore, be arbitrarily small and apparently the light slepton-sneutrino scenario can be accommodated.

Both the above solutions, however, lead to an unstable electroweak symmetry breaking vacuum as the potential becomes unbounded from below. In fact in the case of the AMSB model the requirement of vacuum stability leads to a lower bound on the slepton mass which after allowing for all theoretical uncertainties is approximately 300 GeV. Of course one can argue that we are living in a false vacuum with a life time larger than the age of the universe which makes the requirement of vacuum stability redundant. However, it is difficult to accept this solution uncritically. First of all the calculation of the probability of tunnelling from the false vacuum to the true one, which is rather straightforward in a model with a single scalar, becomes far too complicated in the MSSM with multiple scalars. In fact the earlier calculations of this probability have been criticized by more recent ones (see the second paper of Ref.). Yet there is no way of testing the reliability of the recent calculations as they cannot be verified experimentally. Furthermore tunnelling being a probabilistic phenomenon, the unpleasant possibility that charge and colour symmetry will be broken at the very next moment always remains open. In our opinion, therefore, the false vacuum scenario should remain as a theoretical curiosity unless the complete determination of the sparticle spectrum in future experiments points unmistakably to a set of SUSY parameters leading to an unstable electroweak symmetry breaking vacuum.

In view of the above discussions it is prudent to look for a model with light sleptons without jeopardizing the stability of the vacuum. Such a scenario arises when an SO(10) SUSY GUT directly breaks down to the SM gauge group. Certain U(1) symmetry breaking $`D`$-terms at $`M_G`$ can then lead to a sparticle spectrum with right sleptons naturally heavier than the left slepton/sneutrinos. In fact it was shown in Ref. that in this model one may have sneutrinos in the mass range preferred by the precision electroweak data while the right selectron mass is beyond the kinematic reach of LEP.

This is not to suggest that the above model is the only or even the most appealing model with light sneutrinos. It simply demonstrates that the physics at $`M_G`$ or $`M_P`$ involves too many uncertainties to determine precisely from the boundary conditions at $`M_G`$, what the sparticle spectrum at low energies should look like. Any mass spectrum apparently preferred by the data should, therefore, be taken seriously irrespective of being favoured or disfavoured by the currently prevailing theoretical prejudices. However, further direct experimental tests to confirm the spectrum is urgently needed. We shall now turn our attention to this task.

If the sneutrino and other sparticle masses are indeed in the range preferred by the precision electroweak data then the $`\stackrel{~}{\nu }`$ would be the next lightest supersymmetric particle (NLSP), the lightest neutralino ($`\stackrel{~}{\chi }_1^0`$) being the lightest supersymmetric particle (LSP). Such a sneutrino will decay into the invisible mode $`\stackrel{~}{\nu }\nu \stackrel{~}{\chi }_1^0`$ with 100% branching ratio (BR). As such sneutrinos will act as carriers of missing energy just like the $`\stackrel{~}{\chi }_1^0`$ in a $`R`$-parity conserving model, they have been called virtual lightest supersymmetric particle (VLSP) or the effective lightest supersymmetric particle (ELSP). Such sneutrino pairs produced in association with an energetic photon will lead to the signal $`e^+e^{}\gamma +`$ missing energy This signal if seen over the SM background may indicate the existence of light sneutrinos. In the MSSM there are, however, other processes leading to the same final state. Of course the processes $`e^+e^{}\gamma \stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0`$ is present in all versions of the $`R`$-parity conserving model. If the second lightest neutralino ($`\stackrel{~}{\chi }_2^0`$) decays into the channel $`\stackrel{~}{\chi }_2^0\stackrel{~}{\nu }\nu `$ with a large BR, then all processes belonging to the class $`e^+e^{}\gamma \stackrel{~}{\chi }_i^0\stackrel{~}{\chi }_j^0`$ (i = 1,2), will also contribute to the signal. The actual value of this BR, however, is more model dependent as it depends also on $`m_{\stackrel{~}{\mathrm{}}_R}`$.

The cross sections of all supersymmetric contributions to $`e^+e^{}\gamma +`$ missing energy were exactly computed in Ref.. It was shown that with special kinematic cuts (see below) the signal can be seen over the SM background. However, the results of were computed on the basis of 500 pb-1 of integrated luminosity collected at $`\sqrt{s}=190`$ GeV. More recently LEP has run over a more varied range of $`\sqrt{s}`$ including energies considerably higher than $`\sqrt{s}=190`$ GeV. The total accumulated luminosity also happens to be much larger than that considered in Ref.. Since the sneutrino VLSP signal assumes new significance in the light of the precision electroweak data, a re-examination of the analysis of Ref. is called for.

It may be recalled that all the LEP groups have extensively studied the single $`\gamma `$ \+ missing energy signal. In most of the analyses the cross section of the SM process $`e^+e^{}\nu \overline{\nu }\gamma `$ was measured and the result was then used for neutrino counting. The kinematical cuts were optimized to suppress the backgrounds to this process (radiative Bhabha scattering etc.). This process receives a large contribution when the photon energy is such that an on shell $`Z`$ is produced which subsequently decays to a $`\nu \overline{\nu }`$ pair. An appropriate upper cut on the photon energy excludes this radiative return to the $`Z`$-pole and drastically reduces this cross section. However, in the presence of initial state radiations the photon energy appropriate for radiative return is somewhat smeared and the efficiency of the cut is reduced to some extent (see below). Nevertheless this cut is very useful to improve the signal to background ratio. In a limited number of new physics searches such cuts were also employed by some of the LEP collaborations. For example DELPHI looked for the production of a pair of superlight gravitinos in association with a photon, where $`E_\gamma `$ was restricted to $`E_\gamma <`$ 50 GeV. The dynamics and kinematics of this process is, however, quite different from the production of a pair of heavy sneutrinos and their conclusions cannot be extended to the case under study in a straightforward way.

In the following we shall repeat the analyses of Ref. for realistic energies and luminosities. We, however, find that the data collected at a particular energy or by a particular group is not enough to produce a signal to background ratio which is statistically significant. On the other hand, if the data from all the groups at different energies are combined more significant results may be obtained. Admittedly this is a difficult task and should be carried out by the experts. In particular combining the systematics of different experiments requires special care. Moreover, the integrated luminosities accumulated by each group at a given energy are not always available. For example, in the energy range $`\sqrt{s}`$ = 191.6 - 201.6 GeV the total luminosity collected by OPAL is said to be 212 pb-1. We could not find out further break-ups. On the other hand the break-ups of the total luminosity collected by the other groups at different energies in this range were available in most cases (see Table 1).

The purpose of our simple minded analysis is, therefore, not to obtain rigorous bounds. We are rather interested in illustrating the sensitivity of the existing data to sneutrino masses in the range preferred by the precision electroweak data. We, therefore, often take recourse to various approximations. For example, we roughly estimate the OPAL luminosity at a particular $`\sqrt{s}`$, by scaling the corresponding quantity given by another group, by the ratio of the total luminosity collected by OPAL and the other group in the entire energy range. Whenever the detailed break-ups were not available, we employed this approximation. The information that we could gather from the literature is summarized in Table 1:

| $`\sqrt{s}`$(GeV) | ALEPH | DELPHI | L3 | OPAL | Total $``$ (pb-1) |
| --- | --- | --- | --- | --- | --- |
| 188.6 | 173.6 | 154.7 | 176.4 | 177.3 | 682.0 |
| 191.6 | 28.9 | 25.9 | 29.5 | $``$ 25.7 | $``$ 110.0 |
| 195.5 | 79.9 | 76.4 | 83.0 | $``$ 73.0 | $``$ 312.3 |
| 199.5 | 87.0 | 83.4 | 82.1 | $``$ 76.8 | $``$ 329.3 |
| 201.6 | 44.4 | 40.6 | 36.9 | $``$ 36.4 | $``$ 158.3 |
| 203.7 | $``$ 6.9 | 8.4 | $``$ 8.5 | $``$ 6.3 | $``$ 30.1 |
| 205.2 | 75.3 | 76.2 | $``$ 77.9 | $``$ 57.8 | $``$ 287.2 |
| 206.7 | 122.6 | 121.6 | $``$ 125.6 | $``$ 93.2 | $``$ 463.0 |
| 208.2 | 9.4 | 8.3 | $``$ 9.1 | $``$ 6.8 | $``$ 33.6 |
| Total | 628.0 | 595.5 | $``$ 629.0 | $``$ 553.3 | 2405.8 |

Table 1. The breakups for luminosities accumulated by the LEP groups in different $`\sqrt{s}`$-bins in the range $`188.6<\sqrt{s}<208.2`$ GeV. The approximated luminosities are indicated explicitly and follow from the treatment explained in the text. In Table 2 we present a sample analysis. For this analysis we assume gaugino mass unification and the slepton masses to be flavour independent. The SUSY parameters are chosen to be $`m_{\stackrel{~}{\nu }}`$ = 56 GeV, $`M_2`$ (the SU(2) gaugino mass) = 110 GeV, $`\mu `$ (the higgsino mass parameter) = $`300`$ GeV and $`\mathrm{tan}\beta `$ (the ratio of the vacuum expectation values of the two Higgs bosons) = 10. These immediately lead to the following mass pattern: $`m_{\stackrel{~}{\mathrm{}}_L}`$ (the mass of the left slepton) = 96.96 GeV, $`m_{\stackrel{~}{\chi }^\pm }`$ (the mass of the lighter chargino) = 106.5 GeV, $`m_{\stackrel{~}{\chi }_1^0}`$ (mass of the lightest neutralino) = 55.1 GeV and $`m_{\stackrel{~}{\chi }_2^0}`$ (the mass of the second lightest neutralino) = 106.3 GeV. The mass of the right slepton $`(m_{\stackrel{~}{\mathrm{}}_R})`$ cannot be fixed without further assumption. We assume it to be 100 GeV. As discussed above $`\stackrel{~}{\chi }_2^0`$ can decay into charged leptons and sleptons and, hence, may not decay into the invisible channel $`\nu \stackrel{~}{\nu }`$ with 100% BR. In order to obtain a conservative estimate we, therefore, ignore the $`\stackrel{~}{\chi }_1^0+\stackrel{~}{\chi }_2^0`$ \+ $`\gamma `$ events, although a significant fraction of such events may also lead to the signal under study. Thus in practice the signal in Table 2 can be somewhat larger than what has been shown.

The kinematic cuts imposed are as follows. Only hard photons emitted into the angular interval $`15^{}<\theta _\gamma <165^{}`$ are considered, where $`\theta _\gamma `$ is the polar angle of the photon with respect to the beam direction. Roughly speaking this covers the barrel and the end cap regions of a typical LEP detector. Of course, the detailed geometries are different for different experiments. The detection efficiency of photons, which is detector dependent, is taken to be unity for simplicity. There is a lower cut on the photon energy $`E_\gamma `$ (or $`p_{T_\gamma }`$) which along with the above angular cuts reduces the background from radiative Bhabha scattering and other processes. This lower cut on photon energy is parametrized by a variable $`x_\gamma `$ defined to be $`x_\gamma =E_\gamma /E_{beam}`$. Typically the LEP experiments require $`x_\gamma >0.06`$ in the barrel regions and $`x_\gamma >0.1`$ in the end cap regions. As we are trying to include the end cap region in our analysis we conservatively set $`x_\gamma >0.09`$ for all $`\theta _\gamma `$ which is likely to compensate for the gaps actually present between the end cap and the barrel regions that we neglect in our analysis. In addition there should be an upper cut on the photon energy to exclude events with radiative return to the $`Z`$-pole. Again, in principle, this cut should be a function of $`\sqrt{s}`$ and $`m_{\stackrel{~}{\nu }}`$ to be probed. But in our simple minded analysis the photon energy is restricted to the range $`E_{min}<E_\gamma <60`$ GeV at all c.m. energies. The signal to background ratio can be further improved by optimizing the upper cut with beam energy and $`m_{\stackrel{~}{\nu }}`$.

| $`\sqrt{s}`$ | $`\sigma _{\nu \overline{\nu }\gamma }`$ | $`\sigma _{\stackrel{~}{\nu }\stackrel{~}{\nu }^{}\gamma }`$ | $`\sigma _{\chi _1^0\chi _1^0\gamma }`$ | $``$ | Number | Number |
| --- | --- | --- | --- | --- | --- | --- |
| (GeV) | (pb) | (pb) | (pb) | (pb-1) | of BG Events | of Signal Events |
| 188.6 | 1.303(2.256) | 0.101 | 0.018 | 682 | 886(1534) | 81 |
| 191.6 | 1.268(2.162) | 0.103 | 0.019 | 110 | 139(238) | 13 |
| 195.5 | 1.251(2.040) | 0.105 | 0.019 | 312 | 390(636) | 39 |
| 199.5 | 1.222(1.949) | 0.106 | 0.020 | 329 | 402(641) | 41 |
| 201.6 | 1.208(1.894) | 0.107 | 0.020 | 158 | 188(295) | 20 |
| 203.7 | 1.201(1.851) | 0.107 | 0.020 | 30 | 36(56) | 4 |
| 205.2 | 1.193(1.828) | 0.108 | 0.020 | 287 | 338(517) | 36 |
| 206.7 | 1.188(1.794) | 0.108 | 0.021 | 463 | 546(825) | 59 |
| 208.2 | 1.179(1.771) | 0.108 | 0.021 | 33 | 38(57) | 4 |
| Total |  |  |  |  | 2963(4599) | 297 |

Table 2. The signal and the background rates for different values of $`\sqrt{s}`$. The kinematic cuts employed are $`15^{}<\theta _\gamma <165^{}`$, $`E_\gamma ^{min}<E_\gamma <60`$ GeV where $`E_\gamma ^{min}=x_\gamma E_b`$. The SUSY parameters are $`\mu =300`$ GeV, $`\mathrm{tan}\beta =10`$ and $`M_2=2M_1=110`$ GeV, $`m_{\stackrel{~}{\nu }}=56`$ GeV , $`m_{\stackrel{~}{e}_L}=96.96`$ GeV and $`m_{\stackrel{~}{e}_R}=100`$ GeV.

A few remarks on the method of calculation are now in order. It was noted in Ref. that the cross section for $`e^+e^{}\nu \overline{\nu }\gamma `$ calculated by us was smaller than the cross sections reported by other groups. Subsequently we found a sign error in one of the interfering amplitudes and our results are now in exact agreement with the other groups. The cross section corrected for initial state radiation as discussed in Ref., agrees nicely with the result obtained by using the package CompHEP. There was also a sign error in one of the interfering amplitudes for the process $`e^+e^{}\stackrel{~}{\nu }\stackrel{~}{\nu }^{}\gamma `$ which underestimated the signal cross section. Since both the signal and the background were underestimated in Ref., the conclusions were roughly correct albeit somewhat fortuitously. After corrections our result for $`e^+e^{}\stackrel{~}{\nu }\stackrel{~}{\nu }^{}\gamma `$ agrees with that from CompHEP. Also, the cross sections of $`\stackrel{~}{\chi }_1^0`$ pair + $`\gamma `$ production is in agreement with results from CompHEP. All the results in Table 2 have been obtained by CompHEP and cross checked against our corrected codes.

In Table 2 the cross sections are given for SUSY parameters and kinematic cuts as discussed above. In column 2 we present the cross section for $`e^+e^{}\nu \overline{\nu }\gamma `$. The first number corresponds to cross section of the bare process and the number in parenthesis to that corrected for initial state radiations. The latter is somewhat larger than the former for reasons discussed above. The SUSY signal cross sections are presented in columns 3 and 4, while the total integrated luminosity at each energy as given by Table 1, is presented in column 5. It was already noted in Ref. that the signal cross sections are by and large unaffected by correcting for initial state radiation. The number of background events without (with) initial state radiation is given in column 6, while the total number of signal events (obtained from the information in column 3, 4 and 5) is given in column 7. It is quite clear that the signal($`S`$) to $`\sqrt{\mathrm{background}(B)}`$ ratio at any particular energy is not statistically very significant. However, adding the numbers in the last two columns of Table 2 as shown in the last row, we find that $`S/\sqrt{B}`$ is 4.4 if initial state radiation is taken into account.

It may be recalled that the range of the sneutrino mass preferred by the precision electroweak data is 55 GeV $`<m_{\stackrel{~}{\nu }}<`$ 80 GeV. On the face of it, therefore, it appears as though only the lower edge of this mass range can be probed. However, there are some opportunities for improvement. First of all LEP data above the $`Z`$-pole may be considered ( the present analysis is restricted to $`\sqrt{s}>`$ 189 GeV). In view of the approximate nature of the present analysis we made no attempt to optimize the cuts. For example, a fixed upper cut on $`E_\gamma `$ has been employed in the present analysis at all energies. Alternatively cuts suitably tailored for probing a particular $`m_{\stackrel{~}{\nu }}`$ at a specific C.M. energy may be employed. Finally the contribution of $`e^+e^{}\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_2^0\gamma `$, where $`\stackrel{~}{\chi }_2^0`$ also decay invisibly, may further enhance the signal.

There are reasons to be optimistic about this analysis. It is encouraging to note that there is already an attempt to combine the single photon + missing energy data, at all energies above the $`Z`$-pole, of the four LEP groups to constrain some new physics scenarios. From the combined recoil mass (M) distribution it even appears that there is a modest excess of events over the SM background for M$`>`$ 160 GeV. It should be relatively straightforward to extend this analysis to the case of a pair of heavy invisible sneutrinos produced in association with the photon. If no signal is found, some sneutrino masses in the range preferred by the precision data will be excluded. Even this is an improvement over the current weak lower bound on $`m_{\stackrel{~}{\nu }}`$ ($`0.5m_Z`$) obtained from the $`Z`$-pole data. Of course the latter bound is model independent, while one under consideration depends on gaugino mass unification and on assumptions about $`m_{\stackrel{~}{\mathrm{}}_R}`$. Yet for the first time one will probe $`m_{\stackrel{~}{\nu }}`$ through direct production of invisible sneutrino pairs even if sleptons are beyond the kinematic reach of LEP. In that sense the model dependent assumptions employed are quite different from the assumptions implicit in mSUGRA. The bounds obtained by this method will therefore be complimentary to that obtained, e.g., in Ref.. How a completely model independent bound can be obtained by this method will be discussed below.

Since $`m_{\stackrel{~}{\mathrm{}}_R}`$ was somewhat arbitrarily chosen in our analysis, some discussion of the sensitivity of the signal to this parameter is in order. This parameter does not affect the cross section of $`e^+e^{}\stackrel{~}{\nu }\stackrel{~}{\nu }^{}\gamma `$ which is by far the dominant contribution to the signal for low and moderate $`m_{\stackrel{~}{\nu }}`$. However, if $`m_{\stackrel{~}{\nu }}`$ is close to the upper end of the range preferred by the electroweak precision data, the observability of the signal can indeed be affected by $`m_{\stackrel{~}{\mathrm{}}_R}`$ as can be seen from Figure 1 drawn for $`\sqrt{s}`$ = 205 GeV. The upper (lower) curve corresponds to $`m_{\stackrel{~}{\mathrm{}}_R}=m_{\stackrel{~}{\mathrm{}}_L}`$ ($`m_{\stackrel{~}{\mathrm{}}_R}=1.5m_{\stackrel{~}{\mathrm{}}_L}`$). The other SUSY parameters are as in Table 2. The background cross section at this energy is given in Table 2.

One can use the absence of signal to put a lower bound on the mass of any invisible sneutrino in a model independent way. If we do not make any additional assumption like gaugino mass unification then the neutralino contributions to the $`\gamma +`$ missing energy signal is uncorrelated to the sneutrino contribution. The most conservative lower bound on the sneutrino mass can, therefore, be obtained by comparing the data with the sneutrino contribution alone, which is independent of $`m_{\stackrel{~}{\mathrm{}}_R}`$.

For $`\stackrel{~}{\nu }_\mu `$ and $`\stackrel{~}{\nu }_\tau `$ only the $`Z`$-mediated $`s`$-channel diagrams contribute. Hence, the cross sections are model independent and depend on $`m_{\stackrel{~}{\nu }}`$ alone. The $`\stackrel{~}{\nu }_e`$ contribution depends on $`m_{\stackrel{~}{\nu }}`$, $`M_2`$, $`\mu `$ and tan $`\beta `$. The last three parameters determine the chargino mass and its coupling with the sneutrino. The amplitude has two pieces: the $`s`$-channel $`Z`$-exchange terms and the $`t`$-channel chargino exchange terms. The key observation is that the interference among the $`s`$-channel and $`t`$-channel diagrams are always constructive for all choices of $`M_2`$, $`\mu `$ and tan $`\beta `$. Thus the smallest possible signal cross section can be obtained by neglecting the contribution of the $`t`$-channel diagrams. This corresponds to infinite chargino mass or Higgsino dominated charginos, completely decoupled from the sneutrinos. Now one obtains the minimum signal as a function of $`m_{\stackrel{~}{\nu }}`$ only by assuming three generations of degenerate sneutrinos.

Considering only the $`s`$-channel contribution to the signal we find 96(64) events for $`m_{\stackrel{~}{\nu }}`$ = 50(55) GeV for the entire energy range shown in Table 2. Such a small signal will be swamped by the background. Obviously the current data alone is not adequate to significantly strengthen the existing model independent lower bound on $`m_{\stackrel{~}{\nu }}`$. However, the size of the signal at the Next Linear Collider as given in  suggests that if no signal is found a model independent bound may emerge. The inclusion of other SUSY contributions in specific models, can only strengthen this bound.

Light sneutrinos in the mass range indicated by the precision electroweak data may dramatically influence the SUSY search strategies at the upgraded Tevatron. It has been already been noted since long back, that the hadronically quiet trilepton + $`\overline{)}E_T`$ event resulting from $`\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_2^0`$ pairs is the best channel for SUSY search at this machine. In the presence of light sneutrinos, however, many of the $`\stackrel{~}{\chi }_2^0`$s will decay into the invisible channel $`\nu \stackrel{~}{\nu }^{}`$. Thus the trilepton channel may be considerablly weakened or wiped out depending on the BR of the invisible channel, which is model dependent. On the other hand $`\stackrel{~}{\chi }^\pm `$ pair also has a healthy production cross section at the upgraded Tevatron . Since these charginos will decay into the two body channel $`l+\stackrel{~}{\nu }`$ with a large BR, the opposite sign dilepton + $`\overline{)}E_T`$ signal will be enhanced . Suitable kinematical cuts can show this signal over WW and other relevant backgrounds. From the results of  it appears that the entire range of $`m_{\stackrel{~}{\nu }}`$ preferred by the elctroweak data can be probed.

Acknowledgements : AD thanks Professor H. Dreiner and other members of the Theoretical Physics Group, University of Bonn where part of the work was done. He thanks Professor M. Kobel for discussions. Aseshkrishna Datta is supported by French MNERT fellowship and subsequently by CNRS France. He also likes to thank the Theoretical Physics Group of Abdus Salam ICTP, Trieste for a visit during summer 2001 where this project was started.