text:
# Ultrafast spin dynamics and critical behavior in half-metallic ferromagnet : Sr₂FeMoO₆

## Abstract

Ultrafast spin dynamics in ferromagnetic half-metallic compound S$`\mathrm{r}_2\mathrm{FeMoO}_6`$ is investigated by pump-probe measurements of magneto-optical Kerr effect. Half-metallic nature of this material gives rise to anomalous thermal insulation between spins and electrons, and allows us to pursue the spin dynamics from a few to several hundred picoseconds after the optical excitation. The optically detected magnetization dynamics clearly shows the crossover from microscopic photo-induced demagnetization to macroscopic critical behavior with universal power law divergence of relaxation time for wide dynamical critical region.

Control and manipulation of spins by ultrafast optical excitation, which gives rise to photo-induced magnetization change and magnetic phase transitions in dilute magnetic semiconductor quantum structures , doped semiconductors and ferromagnetic metals , have attracted considerable attention. Recent study on the magnetization dynamics in the photo-excited Ni films with nonlinear optical techniques has revealed ultrafast spin process within 50 fs . Strongly correlated electron systems with half-metallic nature, which have perfectly spin polarized conducting electrons at the ground state, are promising candidates for the study of the photo-induced spin dynamics. These materials have been found to possess exotic physical properties such as colossal magnetoresistance, which have strong application potential . The strong coupling between spin, charge and lattice degrees of freedom in strongly correlated systems makes it possible to manipulate the magnetic properties via cooperative effects induced by optical excitation. In particular, the evidence of photo-induced phase transition accompanied with magnetization changes have been recently reported . In order to understand the nature of these phenomena, it is crucial to investigate the temporal evolution of the spin system in the picosecond time scale. Although some attempts have been made by employing pump-probe spectroscopy , to the best of our knowledge direct investigation of the ultrafast spin dynamics in half-metallic materials has not been reported so far. Such an investigation can be carried out by exploiting the time resolved magneto-optical Kerr effect (MOKE), which has been shown to be a powerful tool to study the ultrafast dynamics of magnetization .

In the present paper we report on the ultrafast pump-probe MOKE and reflectivity study of dynamics of spin and electron systems in the ordered double perovskite $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$. Since Fe$`{}_{}{}^{3+}(3d^5;t_{2g}^3e_g^2,S=5/2)`$ and Mo$`{}_{}{}^{5+}(4d^1;t_{2g}^1,S=1/2)`$ couple antiferromagentically via interatomic exchange interaction and the down-spin electron of Mo5+ is considered itinerant (upper panel of Fig. 1), a conducting ferrimagnetic ground state with half-metallic nature is expected for this material. The density functional calculation also shows that the occupied up-spin band mainly consists of Fe $`3d`$ electrons, while the Fermi level exists within the down-spin band composed of Fe $`t_{2g}`$ and Mo $`t_{2g}`$ electrons. The temperature dependence of magnetization and resistivity measurements under the magnetic field have shown the ferromagnetic phase transition with the Curie temperature $`T_C`$ $``$ 410-450 K Also, this material in the form of polycrystalline ceramics shows inter-grain tunneling type giant magnetoresistance at room temperature because of its spin polarization. Another important property of $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ is the enhancement of MOKE due to spin-orbit coupling of $`t_{2g}`$ electrons in the heavy Mo-atom. The strong MOKE signal enables us to investigate the spin dynamics over a wide temporal range from sub-picosecond to nanosecond and to obtain the critical exponent of the relaxation time at the magnetic phase transition.

The MOKE measurements are carried out on single crystal $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$, grown by floating-zone method, in polar Kerr configuration under the magnetic field of 2000 Oe, where the magnetization is nearly saturated at room temperature , utilizing polarization modulation by a piezo-elastic modulator (CaF2). Figure 1 (a) shows the spectral profiles of ellipticity $`\eta `$ and rotation angle $`\theta `$ at room temperature. A very large MOKE signal, one order of magnitude larger than that of the doped manganites is observed. The MOKE signal is proportional to $`fM`$, where $`M`$ is the magnetization and $`f`$ is determined by the complex refractive index at the probe frequency. Correspondingly, the magneto-optical spectra show resonance, known as the plasma enhancement effect , around 1eV, which is close to the plasma edge. The temperature dependence of the $`\eta `$, probed at a photon energy of 0.95 eV, clearly shows the magnetic phase transition at $`450\mathrm{K}`$, which is the $`T_C`$ of the present sample. Since the reflectivity is almost temperature independent in this temperature range, the sample magnetization can be monitored with the $`\eta `$.

For the pump-probe measurements, a Ti:Sapphire regenerative amplifier system (1 kHz repetition rate) with an Optical Parametric Amplifier (OPA) is used as light source. The second harmonic of the amplified pulses with a pulse duration of 200 fs, photon energy of 3.1 eV and a maximum fluence of 90 $`\mu \mathrm{J}/`$cm2 is used as pump pulse, and its energy is close to the charge transfer excitation from O $`2p`$ to Fe/Mo $`t_{2g}`$ band with down-spin (see the upper panel of Fig.1). The probe pulses from the OPA are tuned to 0.95 eV, at which the ellipticity dominates the MOKE signal rather than the rotation effect (Fig.1 (b)). The pump-probe MOKE measurements are also carried out in polar Kerr configuration and the polarization change of the reflected light from the sample is measured by a balanced detection scheme shown in Fig. 2 (a). By synchronizing the chopper for the pump beam with the regenerative amplifier, the balanced signal is detected and analyzed shot by shot using boxcar integrator and A/D converter. The photo-induced Kerr ellipticity change is measured as the difference between magnetization reversal signals, i.e., $`\mathrm{\Delta }\eta _{Kerr}=\frac{1}{2}[\mathrm{\Delta }\eta (M)\mathrm{\Delta }\eta (M)]`$ by changing the sign of the magnetic field in order to eliminate the contribution from the pump induced optical anisotropy. A sensitivity of 10-3 deg is achieved in our measurement system. The signal is observed to be proportional to the pump beam intensity in all pump-probe measurements.

The inset in Fig. 2 (a) shows the transient reflection change $`\mathrm{\Delta }R/R`$, measured at 300 K. It shows a sharp reduction in the reflectivity during the pump pulse duration, followed by a fast relaxation within $`23\mathrm{p}\mathrm{s}`$ (region (1)) and a fairly long time plateau up to few tens of nanosecond (region (2)). The reflectivity returns to the initial state in 1 ms by heat diffusion (region (3)). Figures 2 (b) and (c) show the $`\mathrm{\Delta }R/R`$ and ellipticity change $`\mathrm{\Delta }\eta _{Kerr}`$ for different temperatures, indicating that the temperature dependence is negligible in $`\mathrm{\Delta }R/R`$, while that is significant in $`\mathrm{\Delta }\eta _{Kerr}`$. The temporal evolution of $`\mathrm{\Delta }\eta _{Kerr}`$ up to 500 ps is shown in Fig.3 (a) for different temperatures. One can observe from Fig. 3(a), that below the Curie point, it can be fitted by $`\mathrm{\Delta }\eta _{Kerr}\left(t\right)=\mathrm{\Delta }\eta _{step}+(\mathrm{\Delta }\eta _{max}\mathrm{\Delta }\eta _{step})(1\mathrm{exp}\left[t/\tau _{spin}\right])`$, where $`\mathrm{\Delta }\eta _{step}`$ describes the instantaneous decrease in the $`\mathrm{\Delta }\eta _{Kerr}`$, while $`\mathrm{\Delta }\eta _{max}`$ is the asymptotic value at the quasi-equilibrium state. The signal, $`\mathrm{\Delta }\eta _{max}`$, increases drastically close to $`T_C`$ as shown in Fig.3 (b). Our measurements reveal nearly linear increase of $`\mathrm{\Delta }\eta _{max}`$ with the pump intensity.

The most striking feature is the very slow spin thermalization observed in $`\mathrm{\Delta }\eta _{Kerr}`$ signal in comparison with the electron thermalization observed in transient $`\mathrm{\Delta }R/R`$ data. In $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$, the behavior of electrons is similar to that of ferromagnetic nickel . Specifically, the electron temperature rises rapidly by the optical excitation and it relaxes within 2 to 3 picoseconds to reach quasi-equilibrium temperature, which is 8-10 K, obtained from $`\mathrm{\Delta }R/R`$, higher than the initial temperature. The fast decay of transient reflectivity indicates that the local heat transfer from electron to the lattice system is completed within a few picoseconds, accompanied by the lattice heat-up to reach quasi-equilibrium temperature. On the other hand, the behavior of the spin system in $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ looks very different from the behavior of the electronic system. Specifically, the very slow spin thermalization (see Fig. 2(c)), which is pronounced at higher temperature, indicates the anomalously small heat exchange between electrons and spins in $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$. Such an electron-spin thermal insulation can be attributed to the nature of the half-metallic electronic structure where conducting electrons are perfectly spin polarized in the down-spin band and isolated from the insulating up-spin band as shown schematically in the upper panel of Fig. 1. Thermal motion of electrons around the Fermi level in the spin polarized conduction band does not increase the spin temperature.

From the observed results, we have the following scenario for the temporal evolution of electron, lattice and spin system in $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$. Initially, during the photo-excitation ($`1\mathrm{ps}`$), the electron system is heated-up and rapidly thermalized due to electron-electron interaction. In this first stage the ellipticity shows a sharp decrease ($`\mathrm{\Delta }\eta _{step}`$). In the next stage, the electron system relaxes by its energy transfer to the lattice system. The electron and lattice systems reach quasi-equilibrium state ($`5\mathrm{p}\mathrm{s}`$) by the electron-phonon interaction, leaving the spin system at its initial temperature. After that, the spin slowly relaxes toward this quasi-equilibrium state through weak heat exchange with the reservoir at quasi-equilibrium temperature. Finally, the system returns to the initial state by heat diffusion.

The sharp decrease in the ellipticity is the major feature of the initial stage of the optical relaxation in $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ (see inset of Fig. 2). The ratio $`\mathrm{\Delta }\eta _{step}/\eta `$ shows a weak temperature dependence (see Fig. 3(c)). Since MOKE signal is proportional to $`fM`$, both photo-induced change in the refractive index ($`\mathrm{\Delta }f/f\mathrm{\Delta }R/R)`$ as well as the photo-induced magnetization change $`\left(\mathrm{\Delta }M/M\right)`$ contribute to $`\mathrm{\Delta }\eta _{step}/\eta `$. Though the relative instantaneous changes in the ellipticity and reflectivity are of the same order, $`\mathrm{\Delta }\eta _{step}/\eta \mathrm{\Delta }R/R0.01`$, the subsequent temporal evolution of $`\mathrm{\Delta }\eta _{Kerr}`$ in the picosecond time scale is very different from that of $`\mathrm{\Delta }R`$. This indicates the direct demagnetization by resonant optical excitation. However, as it has been discussed in recent papers on Ni , it is premature to directly connect the MOKE signal with demagnetization in such ultrafast time scale.

We now discuss the temporal evolution of the MOKE signal in the second stage, when the electron and lattice system have reached the quasi-equilibrium (plateau region (2) in the inset of Fig. 2). The dramatic increase in $`\mathrm{\Delta }\eta _{max}`$(see Fig. 3(b)) and the relaxation time $`\tau _{spin}`$ as the temperature approaches $`T_C`$ indicate that the time resolved signal directly reflects the critical behavior of magnetization at the ferromagnetic phase transition. The spin temperature at the quasi-equilibrium, which can be estimated from Fig. 3(b) and the temperature dependence of the $`\eta `$ (Fig. 1(a)), is in good agreement with the electron temperature estimated from $`\mathrm{\Delta }R/R`$.

It is necessary to emphasize that the time-resolved MOKE measurements give us an unique opportunity to study the critical dynamics of spin system independently from other degrees of freedom and obtain the critical characteristics of the ferromagnetic phase transition. The dynamics of the second order phase transition can be described by the dynamical scaling theory , which allows us to relate the critical behavior of the kinetic parameters (e.g. the relaxation time) to the critical exponents of the static parameters (e.g. correlation length) on the both sides of the critical point. The theory predicts that in the vicinity of $`T_C`$, the relaxation time of the order parameter can be described as $`\tau |TT_C|^{z\nu }`$, where $`\nu `$ and $`z`$ denote the critical exponent of the correlation length and the dynamical critical exponent respectively. Figure 4 shows the temperature dependence of $`\tau _{spin}`$ as a function of $`|T^{}T_C|`$, where $`T^{}`$denotes the quasi-equilibrium temperature at the plateau region (temperature region (2) in the inset of Fig.2(a)). One can clearly observe the power law divergence with $`z\nu =1.22\pm 0.06`$ for the spin relaxation time in the vicinity of $`T_C`$. It should be emphasized that the power law behavior is established in the time scale of few tens of picosecond, while the width of the dynamical critical region is much higher than for conventional metals . The theoretical calculation for the three dimensional Ising and Heisenberg models give $`z\nu 1.30`$ and $`1.37`$, respectively, while the two-dimensional Ising model gives $`z\nu 2.165`$. Therefore, our measurements clearly indicate three dimensionality of the spin system in $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$.

We have presented ultrafast spin dynamics in the ordered double perovskite $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ by using the time-resolved MOKE technique. We have observed, for the first time, extremely slow relaxation of spins, thermally insulated from electron and lattice systems due to the half-metal nature of this material. The thermal insulation of spin system provides us a unique opportunity to examine the non-equilibrium spin dynamics near the critical point in a time scale from picosecond to nanosecond range. Crossover from ultrafast microscopic spin relaxation to macroscopic critical behavior has been clearly demonstrated. In the vicinity of the critical point, the spin relaxation time increases as $`|TT_C|^{(1.22\pm 0.06)}`$ , which is consistent with the theoretical prediction for the 3D ferromagnetic system. We also observe very fast decrease in the MOKE signal $`\mathrm{\Delta }\eta _{Kerr}`$ caused by charge transfer optical excitation. Although the underlying physical mechanism of such an ultrafast phenomenon is not established yet, the distinct difference in the temporal profiles of the photo-induced reflectivity and ellipticity suggests the contribution from ultrafast spin dynamics to nonlinear MOKE signal.

The authors are grateful to S. Miyashita, N. Ito, N. Nagaosa, M. Ueda, Yu. P. Svirko and C. Ramkumar for illuminating discussions. This work is supported in part by a grant-in-aid for COE Research from the Ministry of Education, Science, Sports and Culture of Japan and the New Energy and Industrial Technology Development Organization (NEDO).

Figure captions.

Fig. 1. Magneto-optical Kerr measurements on $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ under the magnetic field of 2000 Oe.(a) Kerr rotation (solid line) and ellipticity (dashed line) spectra at 300K. (b) Temperature profile of linear Kerr ellipticity probed at 0.95eV. The upper panel shows the spin configuration of Fe and Mo ions and a schematic of the electronic band structure of $`\mathrm{Sr}_2\mathrm{FeMoO}_6`$ based on the density-functional calculations by Sawada and Terakura (Ref.).

Fig. 2. (a) Schematic of the experimental setup for the pump-probe magneto-optical Kerr measurements. Inset shows the temporal evolution of transient reflection change $`\mathrm{\Delta }R/R`$ from subpicoseconds up to millisecond. Temporal evolution of the transient reflection $`\mathrm{\Delta }R/R`$ (b) and Kerr ellipticity change $`\mathrm{\Delta }\eta _{Kerr}`$ (c) up to 50ps, measured at 200K, 300K and 400K.

Fig. 3. (a) Temporal evolution of photo-induced Kerr ellipticity $`\mathrm{\Delta }\eta _{Kerr}`$ up to 500 ps, measured at various temperatures. Solid lines are the exponential fit. Temperature profiles of $`\mathrm{\Delta }\eta _{max}`$ (b) and rapid component $`\mathrm{\Delta }\eta _{step}`$ normalized to linear Kerr ellipticity $`\eta `$ (c).

Fig. 4. Temperature dependence of spin relaxation time $`\tau _{spin}`$ as a function of $`|T^{}T_C|`$. The solid line is a power law fit $`|T^{}/T_C1|^{z\nu }`$ for the points near $`T_C.`$ The fit returns $`z\nu =1.22\pm 0.06`$.


format: text

text:
# Three-Dimensional Ising Model and Transfer Matrices

## 1 Introduction

Although over a half-century has passed, solving the 3D Ising model exactly is still an open problem. Anyone who claims to solve this model exactly should, at least, evaluate its partition function, internal energy per site, critical temperature and the critical exponents $`\alpha `$ and $`\beta `$ calculated from the relevant heat capacity and magnetization per site individually. In addition, a crucial test, similar to one L.Onsager  did in 1944, to check whether the results are right or wrong, is that one should compare the series expansion coefficients of the internal energy per site at the high temperature limit with the series obtained from some other methods  e.g., the computer graphic method, at least up to the first three or four nonvanishing terms      .

Among the many various methods for deriving the partition function of the 2D Ising model, transfer matrix method is the oldest and original method. However, the generalization of this method to the 3D case has had relatively little discussion. In this paper, we have no ambition to solve this 3D Ising model satisfying all of the requirements mentioned above. Instead, B.Kaufman’s approach  in the 3D Ising model is carried out step by step. Any approximation is avoided if we possibly can. In the following, it is shown that, when a transfer matrix formalism is set up, a spinor representation can work. 2-dimensional rotations in the direct product space and the feature that all of the eigenvalue equations reduce miraculously to only one equation also appear in the 3D Ising model. Even though the final high-temperature expansion series of internal energy per site is not exactly the same as the computer graphic method’s, these two series do have the same structures. This discrepancy may be related to a dilemma between the choice of the directions of the 2-dimensional rotations in order to find the largest eigenvalue and losing the fascinating feature that all of the eigenvalue equations reduce to only one equation. Be it ever not so perfect, we hope this generalization may lay the foundations for further study.

## 2 Transfer Matrices

Let us consider a simple cubic lattice with $`l`$ layers, each has $`m`$ rows and $`n`$ sites per row. So there are N points on the lattice, N=$`mnl`$. Periodic boundary conditions are used. To each lattice point, with integral coordinates $`\tau `$, $`\rho `$, $`\zeta `$, we assign a spin variable s($`\tau `$,$`\rho `$,$`\zeta `$) which takes two values $`\pm `$ 1. The energy of the configuration is given by

$$E(s)=J\underset{\tau =1}{\overset{n}{}}\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}\{s(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )+s(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )+s(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)\}.$$
(1)
J($`>`$0) is the coupling of a pair of neighboring spins. The partition function

$`𝒵`$ $`=`$ $`{\displaystyle \underset{(s)}{}}e^{E(s)/T}`$ (2)
$`=`$ $`{\displaystyle \underset{(s)}{}}{\displaystyle \underset{\tau =1}{\overset{n}{}}}{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{Ks(\tau ,\rho ,\zeta ,)s(\tau +1,\rho ,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)},`$ (3)
is taken over all the $`2^N`$ possible configurations. Here $`KJ/T`$. Now we factor the partition function into terms each involving only two neighboring spins, giving

$`𝒵`$ $`=`$ $`{\displaystyle \underset{s(1,,)}{}}\mathrm{}{\displaystyle \underset{s(n,,)}{}}<s(1,,)𝒱s(2,,)><s(2,,)𝒱s(3,,)>\mathrm{}`$
$`<s(n,,)𝒱s(1,,)>`$
$`=`$ $`Tr𝒱^n,`$ (5)
where the matrix elements of the transfer matrix $`𝒱`$ are

$$<s(\tau ,,)𝒱s(\tau +1,,)>=\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}e^{Ks(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)}.$$
(6)
$`𝒱`$ can be put into a more convenient form by factoring it into the product of simpler matrices,

$`𝒱`$ $`=`$ $`𝒱_3𝒱_2𝒱_1,`$ (7)
$`<s(\tau ,,)𝒱_1s(\tau +1,)>`$ $`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{Ks(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )},`$ (8)
$`<s(\tau ,,)𝒱_2s(\tau +1,)>`$ $`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )}\delta _{s(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )},`$ (9)
$`<s(\tau ,,)𝒱_3s(\tau +1,,)>`$ $`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)}\delta _{s(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}.`$ (10)
The above decomposition may be checked as follows:

$`<s(\tau ,,)𝒱_3𝒱_2𝒱_1s(\tau +1,,)>`$ (13)
$`=`$ $`{\displaystyle \underset{s(\tau +1,,)}{}}{\displaystyle \underset{s^{}(\tau ,,)}{}}{\displaystyle \underset{s^{}(\tau +1,,)}{}}{\displaystyle \underset{s^{^{\prime \prime }}(\tau ,,)}{}}{\displaystyle \underset{s^{^{\prime \prime }}(\tau +1,,)}{}}<s(\tau ,,)𝒱_3s(\tau +1,,)>`$
$`<s(\tau +1,,)s^{}(\tau ,,)><s^{}(\tau ,,)𝒱_2s^{}(\tau +1,,)>`$
$`<s^{}(\tau +1,,)s^{^{\prime \prime }}(\tau ,,)><s^{^{\prime \prime }}(\tau ,,)𝒱_1s^{^{\prime \prime }}(\tau +1,,)>`$
$`<s^{^{\prime \prime }}(\tau +1,,)s(\tau +1,,)>`$
$`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}\{{\displaystyle \underset{s(\tau +1,\rho ,\zeta )}{}}{\displaystyle \underset{s^{}(\tau ,\rho ,\zeta )}{}}{\displaystyle \underset{s^{}(\tau +1,\rho ,\zeta )}{}}{\displaystyle \underset{s^{^{\prime \prime }}(\tau ,\rho ,\zeta )}{}}{\displaystyle \underset{s^{^{\prime \prime }}(\tau +1,\rho ,\zeta )}{}}`$
$`\left(e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)}\delta _{s(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}\right)\left(\delta _{s(\tau +1,\rho ,\zeta )s^{}(\tau ,\rho ,\zeta )}\right)`$
$`\left(e^{Ks^{}(\tau ,\rho ,\zeta )s^{}(\tau ,\rho +1,\zeta )}\delta _{s^{}(\tau ,\rho ,\zeta )s^{}(\tau +1,\rho ,\zeta )}\right)\left(\delta _{s^{}(\tau +1,\rho ,\zeta )s^{^{\prime \prime }}(\tau ,\rho ,\zeta )}\right)`$
$`\left(e^{Ks^{^{\prime \prime }}(\tau ,\rho ,\zeta )s^{^{\prime \prime }}(\tau +1,\rho ,\zeta )}\right)\left(\delta _{s^{^{\prime \prime }}(\tau +1,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}\right)\}`$
$`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{Ks(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho ,\zeta +1)}.`$ (14)
In the above equation, due to periodic boundary conditions, the identity

$$\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}\delta _{s(\tau ,\rho ,\zeta )s^{}(\tau ,\rho ,\zeta )}=\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}\delta _{s(\tau ,\rho +1,\zeta )s^{}(\tau ,\rho +1,\zeta )}$$
(15)
is used. Furthermore, $`𝒱_1,𝒱_2,𝒱_3`$ can be rewritten as a matrix in the direct product space. Observing from (8), let us define a matrix $`a`$ with matrix elements

$$<s(\tau ,\rho ,\zeta )as(\tau +1,\rho ,\zeta )>=e^{Ks(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}.$$
(16)
$`a`$ $`=`$ $`\left(\begin{array}{cc}e^K& e^K\\ e^K& e^K\end{array}\right)`$ (19)
$`=`$ $`e^K\left(I+e^{2K}\sigma _x\right)`$ (20)
$`=`$ $`\left(2sinh2K\right)^{1/2}e^{K^{}\sigma _x}.`$ (21)
$`I`$ is a $`2\times 2`$ unit matrix. $`tanhK^{}=e^{2K},`$ $`tanhK=e^{2K^{}},`$ $`sinh2Ksinh2K^{}=1.`$ To simplify the matrics $`𝒱_1,𝒱_2,𝒱_3`$, we define

$`X_{i,j}`$ $`=`$ $`\underset{\zeta =1}{\underset{}{\stackrel{mtimesofI}{\stackrel{}{\left(II\mathrm{}I\right)}}}}\mathrm{}\underset{\zeta =j}{\underset{}{\left(I\mathrm{}\stackrel{\rho =i}{\stackrel{}{\sigma _x}}\mathrm{}I\right)}}\mathrm{}`$ (22)
$`\underset{\zeta =l}{\underset{}{\left(II\mathrm{}I\right)}}.`$
$`Y_{i,j}`$ and $`Z_{i,j}`$ are also defined similarly by replacing the Pauli matrix $`\sigma _x`$ with $`\sigma _y`$ and $`\sigma _z`$ respectively.

$`𝒱_1`$ $`=`$ $`\stackrel{mltimes}{\stackrel{}{aa\mathrm{}a}}`$ (23)
$`=`$ $`\left(2sinh2K\right)^{ml/2}{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{K^{}X_{\rho ,\zeta }}.`$ (24)
As for $`𝒱_2`$, we introduce another matrix $`b`$ with matrix elements

$$<s(\tau ,\rho ,\zeta )bs(\tau +1,\rho ,\zeta )>=\delta _{s(\tau ,\rho ,\zeta )s(\tau +1,\rho ,\zeta )}e^{Ks(\tau ,\rho ,\zeta )s(\tau ,\rho +1,\zeta )},$$
(25)
$`b`$ $`=`$ $`\left(\begin{array}{cccc}e^K\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & e^K\hfill & 0\hfill & 0\hfill \\ 0\hfill & 0\hfill & e^K\hfill & 0\hfill \\ 0\hfill & 0\hfill & 0\hfill & e^K\hfill \end{array}\right)=e^{K\sigma _z\sigma _z}=e^{K(\sigma _zI)(I\sigma _z)},`$ (30)
$`𝒱_2`$ $`=`$ $`\stackrel{mltimes}{\stackrel{}{bb\mathrm{}b}}`$ (31)
$`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{KZ_{\rho +1,\zeta }Z_{\rho ,\zeta }}.`$ (32)
Similarly, $`𝒱_3`$ is obtained from a matrix c,

$`c`$ $`=`$ $`e^{K(\sigma _z\stackrel{mtimes}{\stackrel{}{I\mathrm{}I}})(\stackrel{mtimes}{\stackrel{}{II\mathrm{}I}}\sigma _z)}.`$ (33)
$`𝒱_3`$ $`=`$ $`\stackrel{mltimes}{\stackrel{}{cc\mathrm{}c}}`$ (34)
$`=`$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}e^{KZ_{\rho ,\zeta }Z_{\rho ,\zeta +1}}.`$ (35)
## 3 Spinor Representation

$`𝒱_1`$, $`𝒱_2`$,$`𝒱_3`$ in $`2^{ml}`$-space can be related to matrices in 2$`ml`$ -spaces via Dirac $`\mathrm{\Gamma }`$ matrices. The process of reducing the dimensions of $`𝒱`$ had been used in the 2D Ising model. Define a set of matrix $`\mathrm{\Gamma }_{\mu ,\zeta }`$ satisfying anticommutation relations

$$\mathrm{\Gamma }_{\mu ,\zeta }\mathrm{\Gamma }_{\mu ^{^{}},\zeta ^{^{}}}+\mathrm{\Gamma }_{\mu ^{^{}},\zeta ^{^{}}}\mathrm{\Gamma }_{\mu ,\zeta }=2\delta _{\mu \mu ^{^{}}}\delta _{\zeta \zeta ^{^{}}},$$
(36)
$$(\mu ,\mu ^{^{}}=1,2\mathrm{}2m;\zeta ,\zeta ^{^{}}=1,2\mathrm{}l).$$

Every $`\mathrm{\Gamma }_{\mu ,\zeta }`$ is a $`2^{ml}`$ $`\times `$ $`2^{ml}`$ matrix. A possible representation of $`\mathrm{\Gamma }_{\mu ,\zeta }`$ is

$`\mathrm{\Gamma }_{1,1}`$ $`=`$ $`Z_{1,1},`$ (37)
$`\mathrm{\Gamma }_{2,1}`$ $`=`$ $`Y_{1,1},`$ (38)
$`\mathrm{\Gamma }_{3,1}`$ $`=`$ $`X_{1,1}Z_{2,1},`$ (39)
$`\mathrm{\Gamma }_{4,1}`$ $`=`$ $`X_{1,1}Y_{2,1},`$ (40)
$`\mathrm{}`$
$`\mathrm{\Gamma }_{2m1,1}`$ $`=`$ $`X_{1,1}X_{2,1}\mathrm{}X_{m1,1}Z_{m,1},`$ (41)
$`\mathrm{\Gamma }_{2m,1}`$ $`=`$ $`X_{1,1}X_{2,1}\mathrm{}X_{m1,1}Y_{m,1},`$ (42)
$`\mathrm{\Gamma }_{1,2}`$ $`=`$ $`(X_{1,1}X_{2,1}\mathrm{}X_{m,1},)Z_{1,2}=U_1Z_{1,2},`$ (43)
$`\mathrm{\Gamma }_{2,2}`$ $`=`$ $`\left(X_{1,1}X_{2,1}\mathrm{}X_{m,1}\right)Y_{1,2}=U_1Y_{1,2},`$ (44)
$`\mathrm{}`$
$`\mathrm{\Gamma }_{1,\zeta }`$ $`=`$ $`U_1U_2\mathrm{}U_{\zeta 1}Z_{1,\zeta },`$ (45)
$`\mathrm{\Gamma }_{2,\zeta }`$ $`=`$ $`U_1U_2\mathrm{}U_{\zeta 1}Y_{1,\zeta },`$ (46)
$`\mathrm{}`$
$`\mathrm{\Gamma }_{2m1,l}`$ $`=`$ $`U_1U_2\mathrm{}U_{l1}Z_{m,l},`$ (47)
$`\mathrm{\Gamma }_{2m,l}`$ $`=`$ $`U_1U_2\mathrm{}U_{l1}Y_{m,l},`$ (48)
where

$$U_\zeta \left(I\mathrm{}I\right)\mathrm{}\left(\underset{\zeta ^{th}block}{\underset{}{\sigma _x\sigma _x\mathrm{}\sigma _x}}\right)\left(I\mathrm{}I\right)\mathrm{}\left(I\mathrm{}I\right).$$
(49)
The total number of the $`\mathrm{\Gamma }`$ matrix is 2$`ml`$. A special $`2^{ml}\times 2^{ml}`$ matrix U is defined as

$`U`$ $`=`$ $`\underset{mltimes}{\underset{}{\sigma _x\sigma _x\mathrm{}\sigma _x}}`$ (50)
$`=`$ $`i^{ml}{\displaystyle \underset{\zeta =1}{\overset{l}{}}}\left({\displaystyle \underset{\mu =1}{\overset{2m}{}}}\mathrm{\Gamma }_{\mu ,\zeta }\right).`$ (51)
U and $`U_\zeta `$ have the following relations:

$$U^2=I,U_\zeta ^2=I,$$
(52)
$$U\left(I+U\right)=I+U,U_\zeta \left(I+U_\zeta \right)=I+U_\zeta ,$$
(53)
$$U\left(IU\right)=UI,U_\zeta \left(IU_\zeta \right)=U_\zeta I,$$
(54)
$$\{U,\mathrm{\Gamma }_{\mu ,\zeta }\}=0,\{U_\zeta ,\mathrm{\Gamma }_{\mu ,\zeta }\}=0,$$
(55)
$$[U,\mathrm{\Gamma }_{\mu ,\zeta }\mathrm{\Gamma }_{\mu ^{^{}},\zeta ^{^{}}}]=0,[U_\zeta ,\mathrm{\Gamma }_{\mu ,\zeta }\mathrm{\Gamma }_{\mu ^{^{}},\zeta }]=0.$$
(56)
By definition of $`\mathrm{\Gamma }`$, we notice that

$$\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta }=Y_{\rho ,\zeta }Z_{\rho ,\zeta }=iX_{\rho ,\zeta },$$
(57)
then

$$𝒱_1=\left(2sinh2K\right)^{ml/2}\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}e^{iK^{}\mathrm{\Gamma }_{2\rho 1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }}.$$
(58)
Simlarly,

$$\mathrm{\Gamma }_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }=iZ_{\rho ,\zeta }Z_{\rho +1,\zeta },$$
(59)
then we have

$`𝒱_2`$ $`=`$ $`{\displaystyle \underset{\zeta =1}{\overset{l}{}}}\{e^{KZ_{m,\zeta }Z_{1,\zeta }}{\displaystyle \underset{\rho =1}{\overset{m1}{}}}e^{KZ_{\rho ,\zeta }Z_{\rho +1,\zeta }}\}`$ (60)
$`=`$ $`{\displaystyle \underset{\zeta =1}{\overset{l}{}}}\{e^{iKU_\zeta \mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }}{\displaystyle \underset{\rho =1}{\overset{m1}{}}}e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho +1,\zeta }}\}.`$ (61)
With the identity,

$$\left(U_\zeta \mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }\right)^2=U_\zeta \mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }U_\zeta \mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }=I,$$
(62)
$`𝒱_2`$ can be rewritten as

$$𝒱_2=\underset{\zeta =1}{\overset{l}{}}\{\left[\frac{1}{2}(I+U_\zeta )e^{iK\mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }}+\frac{1}{2}(IU_\zeta )e^{iK\mathrm{\Gamma }_{1,\zeta }\mathrm{\Gamma }_{2m,\zeta }}\right]\underset{\rho =1}{\overset{m1}{}}e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho +1,\zeta }}\}.$$
(63)
Since $`U_\zeta `$ commutes with $`\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho +1,\zeta }`$, the projection operators $`\frac{1}{2}(I\pm U_\zeta )`$ project $`𝒱_2`$ into $`2^l`$ pieces. For

$$𝒱_3=\underset{\rho =1}{\overset{m}{}}\left(e^{KZ_{\rho ,l}Z_{\rho ,1}}\underset{\zeta =1}{\overset{l1}{}}e^{KZ_{\rho ,\zeta }Z_{\rho ,\zeta +1}}\right),$$
(64)
the situation seems more complicated.

$`Z_{\rho ,\zeta }Z_{\rho ,\zeta +1}`$ $`=`$ $`\left(Z_{\rho ,\zeta }Z_{\rho +1,\zeta }\right)\left(Z_{\rho +1,\zeta }Z_{\rho +2,\zeta }\right)\mathrm{}\left(Z_{m,\zeta }Z_{1,\zeta +1}\right)\left(Z_{1,\zeta +1}Z_{2,\zeta +1}\right)`$ (66)
$`\mathrm{}\left(Z_{\rho 1,\zeta +1}Z_{\rho ,\zeta +1}\right),`$
$`=`$ $`\left(i\mathrm{\Gamma }_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }\right)\left(i\mathrm{\Gamma }_{2\rho +3,\zeta }\mathrm{\Gamma }_{2\rho +2,\zeta }\right)\mathrm{}`$
$`\left(i\mathrm{\Gamma }_{1,\zeta +1}\mathrm{\Gamma }_{2m,\zeta }\right)\mathrm{}\left(i\mathrm{\Gamma }_{2\rho 1,\zeta +1}\mathrm{\Gamma }_{2\rho 2,\zeta +1}\right),`$
$`=`$ $`i^m\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho +2,\zeta }\mathrm{}\mathrm{\Gamma }_{2m,\zeta }\mathrm{\Gamma }_{1,\zeta +1}\mathrm{}\mathrm{\Gamma }_{2\rho 1,\zeta +1},`$ (67)
$`=`$ $`i^m\mathrm{\Gamma }_{2\rho ,\zeta }\left(\mathrm{\Gamma }_{2\rho +1,\zeta }\mathrm{}\mathrm{\Gamma }_{2\rho 2,\zeta +1}\right)\mathrm{\Gamma }_{2\rho 1,\zeta +1},`$ (68)
$`=`$ $`+i\mathrm{\Gamma }_{2\rho ,\zeta }W_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1},`$ (69)
$`=`$ $`iW_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}.`$ (70)
$`W_{2\rho +1,\zeta }`$ is defined as

$`W_{2\rho +1,\zeta }`$ $`=`$ $`i^{m1}\mathrm{\Gamma }_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho +2,\zeta }\mathrm{}\mathrm{\Gamma }_{2\rho 3,\zeta +1}\mathrm{\Gamma }_{2\rho 2,\zeta +1}`$ (71)
$`=`$ $`I\mathrm{}I\underset{m1timesof\sigma _x}{\underset{}{\stackrel{(\rho +1,\zeta )}{\stackrel{}{\sigma _x}}\mathrm{}\stackrel{(\rho 1,\zeta +1)}{\stackrel{}{\sigma _x}}}}I\mathrm{}I.`$ (72)
$`W_{2\rho +1,\zeta }`$ has the property that it anticommutes with $`\mathrm{\Gamma }_{\mu ,\alpha }`$ inside the region that the integral coordinates $`(\mu ,\alpha )`$ from ($`2\rho +1,\zeta `$) to ( $`2\rho 2,\zeta +1`$ ), whereas it commutes with $`\mathrm{\Gamma }_{\mu ,\alpha }`$ outside of that region.

$`Z_{\rho ,l}Z_{\rho ,1}`$ $`=`$ $`(Z_{\rho ,l}Z_{\rho +1,l})\mathrm{}(Z_{m,l}Z_{1,1})(Z_{1,1}Z_{2,1})\mathrm{}(Z_{\rho 1,1}Z_{\rho ,1})`$ (73)
$`=`$ $`i^mU\mathrm{\Gamma }_{2\rho ,l}\mathrm{\Gamma }_{2\rho +1,l}\mathrm{}\mathrm{\Gamma }_{2m,l}\mathrm{\Gamma }_{1,1}\mathrm{\Gamma }_{2,1}\mathrm{}\mathrm{\Gamma }_{2\rho 1,1}`$ (74)
$`=`$ $`+iUW_{2\rho +1,l}\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}.`$ (75)
$`W_{2\rho +1,l}`$ $`=`$ $`i^{m1}\mathrm{\Gamma }_{2\rho +1,l}\mathrm{\Gamma }_{2\rho +2,l}\mathrm{}\mathrm{\Gamma }_{2m,l}\mathrm{\Gamma }_{1,1}\mathrm{}\mathrm{\Gamma }_{2\rho 2,1}`$ (76)
$`=`$ $`\underset{(1,1)}{\underset{}{\sigma _x}}\mathrm{}\underset{(\rho 1,1)}{\underset{}{\sigma _x}}I\mathrm{}I\underset{(\rho +1,l)}{\underset{}{\sigma _x}}\mathrm{}\underset{(m,l)}{\underset{}{\sigma _x}}.`$ (77)
Then we have

$$𝒱_3=\underset{\rho =1}{\overset{m}{}}\{e^{+iKUW_{2\rho +1,l}\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}\underset{\zeta =1}{\overset{l1}{}}e^{iKW_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}}\}.$$
(78)
A remarkable observation of Kaufman is that decomposing $`𝒱`$ of the 2D Ising model into the product of factors like $`e^{\frac{\theta }{2}\mathrm{\Gamma }\mathrm{\Gamma }}`$, which is interpreted as a two-dimensional rotation with rotation angle $`\theta `$ in the direct product space. We follow this spirit and decompose the factors into several 2D rotations,

$`e^{+iKUW_{2\rho +1,l}\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}`$
$`=`$ $`{\displaystyle \frac{1}{2}}(I+U)e^{+iKW_{2\rho +1,l}\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}+{\displaystyle \frac{1}{2}}(IU)e^{iKW_{2\rho +1,l}\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}`$
$`=`$ $`{\displaystyle \frac{1}{2}}(I+U)\{{\displaystyle \frac{1}{2}}(I+W_{2\rho +1,l})e^{+iK\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}+{\displaystyle \frac{1}{2}}(IW_{2\rho +1,l})e^{iK\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}\}`$
$`+{\displaystyle \frac{1}{2}}(IU)\{{\displaystyle \frac{1}{2}}(I+W_{2\rho +1,l})e^{iK\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}+{\displaystyle \frac{1}{2}}(IW_{2\rho +1,l})e^{+iK\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}\}`$
$`e^{iKW_{2\rho +1,\zeta }\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}}`$
$`=`$ $`{\displaystyle \frac{1}{2}}(I+W_{2\rho +1,\zeta })e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}}+{\displaystyle \frac{1}{2}}(IW_{2\rho +1,\zeta })e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}}.`$ (81)
In the 2D Ising model the $`𝒱`$ are decomposed into 2 pieces. This is not so simple in the 3D Ising model. Due to the projection operators $`\frac{1}{2}`$ (I $`\pm `$ U), $`\frac{1}{2}`$ (I $`\pm `$ $`U_\zeta `$), $`\frac{1}{2}`$(I $`\pm `$ W), $`𝒱`$ has $`2^l\times (2\times 2^{ml})`$ pieces. Only one piece will produce the largest eigenvalue, which dominates the value of partition function.

Since $`W_{2\rho +1,\zeta }`$ may commute or anticommute with $`\mathrm{\Gamma }_{\mu ,\alpha }`$, we have to check the commutation relations between the projection operators and the product of the factors like $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$, since the dogma states that two matrices are simultaneously diagonized if and only if these two matrices commute. $`𝒱`$ comprises the product of projection operators and a lot of $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$. So if the commutation relations are not valid, the whole scheme may break down. Commuting $`\frac{1}{2}(I\pm U)`$ with the product of all factors like $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$ does not give any trouble. Any single projection operator $`\frac{1}{2}(I\pm U_\zeta )`$ or $`\frac{1}{2}(I\pm W)`$ may do not commute with some $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$. However, fortunately, the product of all the projections in any one piece of $`𝒱`$, e.g.

$$P=\frac{1}{2}(I+U)\left(\underset{\zeta ^{}=1}{\overset{l}{}}\frac{1}{2}(I+U_\zeta ^{})\right)\left(\underset{\rho =1}{\overset{m}{}}\underset{\zeta =1}{\overset{l}{}}\frac{1}{2}(I+W_{2\rho 1,\zeta })\right),$$
(82)
do commute with any $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$, because $`e^{\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }}`$ passes through P will change the sign of $`\frac{1}{2}\theta \mathrm{\Gamma }\mathrm{\Gamma }`$ even number times and it does not change sign eventually.

How to choose the proper piece and obtain the largest eigenvalue? We have no answer. A possible rule may be followed that all eigenvalue equations should reduce to one equation. We will show this point in the next section. This is a beautiful feature in the 2D Ising model. Let us consider one possible piece of $`𝒱`$,

$`\stackrel{~}{𝒱}`$ $`=`$ $`(2sinh2K)^{ml/2}\stackrel{~}{𝒱}^{^{}},`$ (83)
$`\stackrel{~}{𝒱}^{^{}}`$ $`=`$ $`\{{\displaystyle \underset{\rho =1}{\overset{m}{}}}e^{+iK\mathrm{\Gamma }_{2\rho 1,1}\mathrm{\Gamma }_{2\rho ,l}}{\displaystyle \underset{\zeta =1}{\overset{l1}{}}}e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho 1,\zeta +1}}\}`$ (84)
$`\{{\displaystyle \underset{\zeta ^{}=1}{\overset{l}{}}}e^{iK\mathrm{\Gamma }_{1,\zeta ^{}}\mathrm{\Gamma }_{2m,\zeta ^{}}}{\displaystyle \underset{\rho ^{}=1}{\overset{m1}{}}}e^{iK\mathrm{\Gamma }_{2\rho ,\zeta }\mathrm{\Gamma }_{2\rho +1,\zeta }}\}`$
$`\{{\displaystyle \underset{\zeta ^{^{\prime \prime }}=1}{\overset{l}{}}}{\displaystyle \underset{\rho ^{^{\prime \prime }}=1}{\overset{m}{}}}e^{ik^{}\mathrm{\Gamma }_{2\rho ^{^{\prime \prime }}1,\zeta ^{^{\prime \prime }}}\mathrm{\Gamma }_{2\rho ^{^{\prime \prime }},\zeta ^{^{\prime \prime }}}}\}.`$
In essence, $`\stackrel{~}{𝒱}`$ includes the repetition $`l`$ times of the same rotations as in the 2D Ising model, appearing in the second and third brackets of (79), and the new rotations in the first bracket of (79), relating to the third dimensional coupling beyond the 2D Ising model.

## 4 Eigenvalue Equations

The rotation operator in the spinor representation,

$$S_{\lambda \sigma }(\theta )=e^{\frac{1}{2}\theta \mathrm{\Gamma }_\lambda \mathrm{\Gamma }_\sigma }(\lambda \sigma ),$$
(85)
has a one-to-one correspondence to the 2D rotational matrix $`\omega (\lambda \sigma \theta )`$ of the $`\mathrm{\Gamma }`$ matrix.

$`S_{\lambda \sigma }^1(\theta )\mathrm{\Gamma }_\alpha S_{\lambda \sigma }(\theta )`$ $`=`$ $`{\displaystyle \underset{\kappa }{}}\omega (\lambda \sigma \theta )_{\alpha \kappa }\mathrm{\Gamma }_\kappa ,`$ (86)
$`S_{\lambda \sigma }^1(\theta )\mathrm{\Gamma }_\lambda S_{\lambda \sigma }(\theta )`$ $`=`$ $`\mathrm{\Gamma }_\lambda \mathrm{cos}\theta +\mathrm{\Gamma }_\sigma \mathrm{sin}\theta ,(\lambda \sigma ),`$ (87)
$`S_{\lambda \sigma }^1(\theta )\mathrm{\Gamma }_\sigma S_{\lambda \sigma }(\theta )`$ $`=`$ $`\mathrm{\Gamma }_\lambda \mathrm{sin}\theta +\mathrm{\Gamma }_\sigma \mathrm{cos}\theta ,(\sigma \lambda ),`$ (88)
$`S_{\lambda \sigma }^1(\theta )\mathrm{\Gamma }_\alpha S_{\lambda \sigma }(\theta )`$ $`=`$ $`\mathrm{\Gamma }_\alpha ,(\alpha \lambda ,\alpha \sigma ).`$ (89)
$`S_{\lambda \sigma }`$($`\theta `$) is a 2D rotations in the direct product space. The rotations, $`e^{\pm \frac{\theta }{2}\mathrm{\Gamma }\mathrm{\Gamma }}`$, have different rotational angles, $`+\theta `$ and $`\theta `$. We mean they have different directions of rotations. The eigenvalues of the rotational matrix $`\omega `$ are 1 with $`(2ml2)`$-fold degeneracies and $`e^{\pm i\theta }`$ two nondegenerate eigenvalues, whereas the eigenvalues of $`S_{\lambda \sigma }`$ in the spinor representation are $`e^{\pm i\frac{\theta }{2}}`$ each with $`2^{ml1}`$ -fold degeneracies.

The correspondence with $`\stackrel{~}{𝒱}`$ is

$`\stackrel{~}{\omega }`$ $`=`$ $`(2sinh2K)^{ml/2}\stackrel{~}{\omega }^{^{}}`$ (90)
$`\stackrel{~}{\omega }^{^{}}`$ $`=`$ $`\{{\displaystyle \underset{\rho =1}{\overset{m}{}}}\omega (2\rho 1,1;2\rho ,l+2iK){\displaystyle \underset{\zeta =1}{\overset{l1}{}}}\omega (2\rho ,\zeta ;2\rho 1,\zeta +12iK)\}`$ (91)
$`\{{\displaystyle \underset{\zeta ^{^{}}=1}{\overset{l}{}}}\omega (1,\zeta ^{^{}};2m,\zeta ^{^{}}2iK){\displaystyle \underset{\rho ^{^{}}=1}{\overset{m1}{}}}\omega (2\rho ,\zeta ;2\rho +1,\zeta 2iK)\}`$
$`\{{\displaystyle \underset{\zeta ^{^{\prime \prime }}=1}{\overset{l}{}}}{\displaystyle \underset{\rho ^{^{\prime \prime }}=1}{\overset{m}{}}}\omega (2\rho ^{^{\prime \prime }}1,\zeta ^{^{\prime \prime }};2\rho ^{^{\prime \prime }},\zeta ^{^{\prime \prime }}2iK^{})\}`$
$`=`$ $`\omega _3\omega _2\omega _1.`$ (92)
$`\omega _1^{\frac{1}{2}}\omega _2\omega _1^{\frac{1}{2}}`$ $`=`$ $`\left(\begin{array}{ccc}\begin{array}{cc}\mathrm{\Omega }& \\ & \mathrm{\Omega }\end{array}\hfill & \begin{array}{cc}& \end{array}\hfill & \hfill \\ \begin{array}{cc}& \end{array}\hfill & \begin{array}{cc}\mathrm{}& \\ & \mathrm{}\end{array}\hfill & \begin{array}{cc}& \end{array}\hfill \\ \hfill & \begin{array}{cc}& \end{array}\hfill & \begin{array}{cc}\mathrm{\Omega }& \\ & \mathrm{\Omega }\end{array}\hfill \end{array}\right)_{2ml\times 2ml},`$ (112)
$`\mathrm{\Omega }`$ $`=`$ $`\left(\begin{array}{c}A\\ B^{}\\ 0\\ \\ 0\\ B\end{array}\begin{array}{c}B\\ A\\ B^{}\\ \\ \mathrm{}\\ 0\end{array}\begin{array}{c}0\\ B\\ A\\ \\ 0\\ \mathrm{}\end{array}\begin{array}{c}\mathrm{}\\ 0\\ B\\ \mathrm{}\\ B^{}\\ 0\end{array}\begin{array}{c}0\\ \mathrm{}\\ 0\\ \\ A\\ B^{}\end{array}\begin{array}{c}B^{}\\ 0\\ \mathrm{}\\ \\ B\\ A\end{array}\right)_{2m\times 2m},`$ (149)
$`A`$ $`=`$ $`\left(\begin{array}{cc}c^{}c& is^{}c\\ is^{}c& c^{}c\end{array}\right),B=\left(\begin{array}{cc}\frac{1}{2}& is(\frac{1+c^{}}{2})\\ is(\frac{1+c^{}}{2})& \frac{1}{2}\end{array}\right),`$ (154)
$`B^{}`$ $`=`$ $`\left(\begin{array}{cc}\frac{1}{2}& is(\frac{1+c^{}}{2})\\ is(\frac{1+c^{}}{2})& \frac{1}{2}\end{array}\right),`$ (157)
$$s\mathrm{sinh}2K,c\mathrm{cosh}2K,s^{}\mathrm{sinh}2K^{},c^{}\mathrm{sinh}2K^{}.$$
(158)
The matrix $`\mathrm{\Omega }`$ is just the same matrix considered in the 2D Ising model. $`\omega _1^{\frac{1}{2}}\omega _2\omega _1^{\frac{1}{2}}`$ is in symmetric form such that its eigenvalue equations are much more easier to handle.

$`\omega _3`$ $``$ $`{\displaystyle \underset{\rho =1}{\overset{m}{}}}\omega (2\rho 1,1;2\rho ,l2iK){\displaystyle \underset{\zeta =1}{\overset{l1}{}}}\omega (2\rho ,\zeta ;2\rho 1,\zeta +12iK)`$ (159)
$`=`$ $`\left(\begin{array}{c}𝒜\\ ^{}\\ 0\\ \\ 0\\ \end{array}\begin{array}{c}\\ 𝒜\\ ^{}\\ \\ \mathrm{}\\ 0\end{array}\begin{array}{c}0\\ \\ 𝒜\\ \\ 0\\ \mathrm{}\end{array}\begin{array}{c}\mathrm{}\\ 0\\ \\ \mathrm{}\\ ^{}\\ 0\end{array}\begin{array}{c}0\\ \mathrm{}\\ 0\\ \\ 𝒜\\ ^{}\end{array}\begin{array}{c}^{}\\ 0\\ \mathrm{}\\ \\ \\ 𝒜\end{array}\right)_{2ml\times 2ml},`$ (196)
$`𝒜`$ $`=`$ $`\left(\begin{array}{c}c\\ \end{array}\begin{array}{c}\\ c\end{array}\begin{array}{c}\\ \\ \mathrm{}\end{array}\begin{array}{c}\\ \\ \\ c\end{array}\right)_{2m\times 2m},^{}=\left(\begin{array}{c}0\\ 0\\ 0\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}is\\ 0\\ 0\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}0\\ 0\\ 0\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ is\\ 0\\ \mathrm{}\\ \end{array}\begin{array}{c}\\ \\ 0\\ 0\\ \\ \end{array}\begin{array}{c}\\ \\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \end{array}\begin{array}{c}\\ \\ \\ \\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \\ \\ 0\\ 0\end{array}\begin{array}{c}\\ \\ \\ \\ \\ \\ is\\ 0\end{array}\right)_{2m\times 2m},`$ (275)
$``$ $`=`$ $`\left(\begin{array}{c}0\\ is\\ 0\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}0\\ 0\\ 0\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}\mathrm{}\\ 0\\ 0\\ is\\ \\ \end{array}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ 0\\ 0\\ \\ \end{array}\begin{array}{c}\\ \\ 0\\ 0\\ \mathrm{}\\ \end{array}\begin{array}{c}\\ \\ \mathrm{}\\ \mathrm{}\\ \\ \end{array}\begin{array}{c}\\ \\ \\ \\ \\ \\ 0\\ is\end{array}\begin{array}{c}\\ \\ \\ \\ \\ \\ 0\\ 0\end{array}\right)_{2m\times 2m}.`$ (334)
Now we proceed to solve the eigenvalue equation

$$\stackrel{~}{\omega }\mathrm{\Psi }=\lambda \mathrm{\Psi },$$
(335)
where $`\stackrel{~}{\omega }=(2sinh2K)^{ml/2}\omega _3\left[\omega _1^{1/2}\omega _2\omega _1^{1/2}\right]`$,

$$\mathrm{\Psi }=\left(\begin{array}{c}z\psi _0\\ z^2\psi _0\\ z^3\psi _0\\ \mathrm{}\\ z^l\psi _0\end{array}\right),\psi _0=\left(\begin{array}{c}yu\\ y^2u\\ y^3u\\ \mathrm{}\\ y^mu\end{array}\right),u=\left(\begin{array}{c}u_1\\ u_2\end{array}\right).$$
(336)
By imposing the constraint,

$$z^l=1,orz=e^{i\frac{\pi t_2}{l}}(t_2=1,3,5\mathrm{}),$$
(337)
one reduces eigenvalue equation (97) to

$$\left(𝒜+z+z^1^{}\right)\mathrm{\Omega }\psi _0=\lambda \psi _0.$$
(338)
Further, imposing the constraint,

$$y^m=1,ory=e^{i\frac{\pi t_1}{m}}(t_1=\mathrm{\hspace{0.33em}1},3,5\mathrm{}),$$
(339)
(100) is reduced to

$$𝒟\left(A+yB+y^1B^{}\right)u=\lambda u,$$
(340)
where

$$𝒟=\left(\begin{array}{cc}c& iz^1s\\ izs& c\end{array}\right),$$
(341)
$`A+yB+y^1B^{}`$ (342)
$`=`$ $`\left(\begin{array}{cc}c^{}c\mathrm{cos}\frac{\pi t_1}{m}& is^{}ci(\frac{1+c^{}}{2})se^{\frac{i\pi t_1}{m}}i(\frac{1+c^{}}{2})se^{i\frac{\pi t_1}{m}}\\ is^{}c+i(\frac{1+c^{}}{2})se^{i\frac{\pi t_1}{m}}+i(\frac{1+c^{}}{2})se^{i\frac{\pi t_1}{m}}& c^{}c\mathrm{cos}\frac{\pi t_1}{m}\end{array}\right).`$ (345)
The eigenvalue equation (102) is further reduced to

$$\lambda ^2\mathrm{\hspace{0.33em}2}\mathrm{cosh}\gamma \lambda +\mathrm{\hspace{0.33em}1}=\mathrm{\hspace{0.33em}0}.$$
(346)
$`\lambda =\mathrm{exp}\left(\pm \gamma \right)`$ is the solution of $`\lambda `$. $`\gamma `$ is determined by

$$\mathrm{cosh}\gamma =\frac{c^3}{s}c(\mathrm{cos}\theta _1+\mathrm{cos}\theta _2)+sc\mathrm{cos}\theta _1\mathrm{cos}\theta _2+s^2\mathrm{sin}\theta _1\mathrm{sin}\theta _2,$$
(347)
where two continuous variables, $`\theta _1`$, $`\theta _2`$, are obtained by taking the thermodynamic limit, $`m,l\mathrm{}`$, $`\frac{\pi t_1}{m}\theta _1`$ , $`\frac{\pi t_2}{l}\theta _2`$.

The partition function of $`𝒱`$ is

$$𝒵\lambda ^N.$$
(348)
The free energy per site under the thermodynamic limit is

$`f`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}\underset{N\mathrm{}}{lim}N^1ln𝒵`$ (349)
$`=`$ $`{\displaystyle \frac{1}{\beta }}\mathrm{ln}(2\mathrm{sinh}2K)^{1/2}{\displaystyle \frac{1}{8\pi ^2\beta }}{\displaystyle _0^{2\pi }}{\displaystyle _0^{2\pi }}\gamma (\theta _1,\theta _2)d\theta _1d\theta _2,`$ (350)
where $`\beta \frac{1}{T}.`$ For convenience, $`J1`$, the internal energy per site is

$`u`$ $`=`$ $`{\displaystyle \frac{}{\beta }}(\beta f)`$ (351)
$`=`$ $`\mathrm{coth}2K{\displaystyle \frac{1}{8\pi ^2}}{\displaystyle _0^{2\pi }}{\displaystyle _0^{2\pi }}{\displaystyle \frac{\gamma }{K}}𝑑\theta _1𝑑\theta _2.`$ (352)
## 5 High Temperature Limit

Let us expand $`u`$ in terms of $`x`$ at high temperature limit, $`x`$ small, and compare the series of $`u`$ obtained from the computer graphic method.

$$x\mathrm{tanh}K,c=\mathrm{cosh}2K=\frac{1+x^2}{1x^2},s=\mathrm{sinh}2K=\frac{2x}{1x^2},$$
(353)
$`\mathrm{cosh}\gamma `$ $`=`$ $`{\displaystyle \frac{\left(1+x^2\right)^3}{\left(1x^2\right)^22x}}{\displaystyle \frac{1+x^2}{1x^2}}(\mathrm{cos}\theta _1+\mathrm{cos}\theta _2)`$ (354)
$`+{\displaystyle \frac{2x\left(1+x^2\right)}{(1x^2)^2}}\mathrm{cos}\theta _1\mathrm{cos}\theta _2+{\displaystyle \frac{4x^2}{\left(1x^2\right)^2}}\mathrm{sin}\theta _1\mathrm{cos}\theta _2,`$
$`{\displaystyle \frac{\gamma }{K}}`$ $`=`$ $`\mathrm{\hspace{0.33em}2}(\mathrm{sinh}\gamma )^1\{{\displaystyle \frac{(1+x^2)^2(110x^2+x^4)}{(1x^2)^24x^2}}{\displaystyle \frac{2x}{1x^2}}(\mathrm{cos}\theta _1+\mathrm{cos}\theta _2)`$
$`+{\displaystyle \frac{1+6x^2+x^4}{(1x^2)^2}}\mathrm{cos}\theta _1\mathrm{cos}\theta _2+{\displaystyle \frac{4x(1+x^2)}{(1x^2)^2}}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\}.`$ (355)
With the help of computer program Maple V Release 5.1, we get our result

$$u=3x8x^328x^5132x^7832x^9+O(x^{11}).$$
(356)
Comparing the result obtained by the computer graphic method, the series expansion of $`u`$ for a simple cubic lattice can be transformed from the partition function,

$$𝒵^{\frac{1}{N}}=2cosh^3K\left(1+\mathrm{\hspace{0.33em}3}x^4+\mathrm{\hspace{0.33em}22}x^6+\mathrm{\hspace{0.33em}187.5}x^8+\mathrm{\hspace{0.33em}1980}x^{10}+O(x^{12})\right),$$
(357)
so

$`u`$ $`=`$ $`{\displaystyle \frac{}{K}}\mathrm{ln}𝒵^{\frac{1}{N}}`$ (358)
$`=`$ $`3x12x^3120x^51332x^717676x^9+O(x^{11}).`$ (359)
(115 ) and (118) have very similar structures. The coefficient of $`x^1`$ does not exist in our result (115) though it exists naively in $`coth2K`$ of (111). The first nonvanishing term $`3x`$ of (115) and the vanishing of all coefficients of even powers of $`x`$ are the same as (118 ). All terms, up to $`O(x^9)`$, with minus signs in (115) are also the same as (118). The reason why (115) and (118) do not have the same first three or four coefficients may come from the fact that we cannot decipher precisely the largest eigenvalue of $`𝒱`$ . On the other hand, if we choose different directions for the 2D rotations, the eigenvalue equations may not be reduced to only one equation. This is a dilemma.

The correct relation of $`\mathrm{cosh}\gamma `$ implies the correct thermodynamic quantities of the 3D Ising model. To guess the correct relation as (106 ) may be a promising way for finding the right resolution of the 3D Ising model. For example, trying to repair (106), if we multiply the matrices $`𝒜`$, $``$ and $`^{}`$, or $`\omega _3`$, with a correcting factor $`\mathrm{\Phi }(x)`$, then we have

$$\mathrm{cosh}\gamma =\frac{c^3}{s}c\mathrm{cos}\theta _1+\mathrm{\Phi }(c\mathrm{cos}\theta _2+sc\mathrm{cos}\theta _1\mathrm{cos}\theta _2+s^2\mathrm{sin}\theta _1\mathrm{sin}\theta _2).$$
(360)
With the help of Maple V, setting

$$\mathrm{\Phi }(x)=1x^2\frac{27}{2}x^4\frac{249}{2}x^6\frac{12325}{8}x^8O(x^{10}),$$
(361)
then the series expansion of $`u`$ has the same result as (118). $`\mathrm{\Phi }(x)`$ may be explained as, something like, a weighting function for different directions of 2D rotations.

## 6 Acknowledgements

One of authors, S.L. Lou, would like to thank H.C. Lee, Friday Lin and C.Y. Lin for many helpful discussions about four years ago.


text:
# Effects of spin-elastic interactions in frustrated Heisenberg antiferromagnets

## Abstract

The Heisenberg antiferromagnet on a compressible triangular lattice in the spin-wave approximation is considered. It is shown that the interaction between quantum fluctuations and elastic degrees of freedom stabilize the low symmetric L-phase with a collinear Néel magnetic ordering. Multi-stability in the dependence of the on-site magnetization is on an uniaxial pressure is found.

There is a long standing interest in properties of frustrated quantum antiferromagnets because they display a large variety of behavior: e.g. conventional spin ordering, spin liquid ground states, or chirality transitions . The antiferromagnetic Heisenberg model on a triangular lattice is the simplest two-dimensional frustrated system. Quasi-two-dimensional triangular lattice antiferromagnets are not rare objects. It was recently shown  that the materials $`AM(SO_4)_2`$ belonging to the Yavapaiite family may be considered as realizations of a quasi-two-dimensional triangular lattice quantum (e.g. M=Ti with S=1/2) and quasiclassical (e.g. M=Fe with S=5/2) antiferromagnets. The rhombohedral $`\beta `$-phase and monoclinic $`\alpha `$-phase of solid oxygen consist of weakly interacting planes (triangular in the case of $`\beta O_2`$ and rectangular in the case of $`\alpha O_2`$ ) whose magnetic properties are described by the easy-plane Heisenberg antiferromagnetic model with spin $`S=1`$ .

Unconventional properties of the $`\kappa (BEDTTTF)_2X`$ family can be described by a Hubbard model on a triangular lattice with an anisotropic interactions . An antiferromagnetic Heisenberg model was also used in  to describe the properties of a triangular monolayer of Pb and Sn adatoms on the (111) surface of Ge form.

From the theoretical point of view the antiferromagnetic Heisenberg model on a triangular lattice has attracted particular attention since Anderson’s proposal of a resonance-valence-bond state for this model . From that time much work has been done to understand the nature of its ground state. There is now an almost consensus based on variety of methods  that the ground state has a conventional three-sublattice order but the situation is close to marginal and quantum fluctuations are important . Quite recently the effect of quantum spin fluctuations on the ground-state properties of spatially anisotropic Heisenberg antiferromagnet in the linear spin-wave approximation was explored . The staggered magnetization and the magnon dispersion were calculated as functions of the ratio of the antiferromagnetic exchange between the second and first neighbors,$`J_2/J_1`$. As a physical tool which allows to tune the ratio $`J_2/J_1`$ Merino et al  invoke an uniaxial stress within a layer whereby the relative distances between atoms are changed.

In this paper we investigate the ground-state properties of the compressible frustrated antiferromagnets, i.e. those systems where the coupling between magnetic and elastic degrees of freedom is taken into account,under the action of uniaxial pressure. We treat the system self-consistently considering the spin subsystem in the linear spin-wave approximation and using the continuum approach for the elastic degrees of freedom. We show that as a result of this method the behavior of the on-site magnetization as a function of the pressure differs significantly from the behavior which was obtained in the framework of the spatially anisotropic Heisenberg model . In particular, the on-site-magnetization does not vanish in the whole interval of stability of the collinear phase.

We consider the two-dimensional triangular lattice. The Hamiltonian of the spin subsystem, interacting with the displacements of the magnetic atoms, is

$`H={\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{n},\stackrel{}{a}}{}}J_\stackrel{}{a}\stackrel{}{S}_\stackrel{}{n}\stackrel{}{S}_{\stackrel{}{n}+\stackrel{}{a}}`$ (1)
Here $`\stackrel{}{S}_\stackrel{}{n}`$ is the spin operator of the atom located in the $`\stackrel{}{n}`$-th site of the triangular lattice: $`\stackrel{}{n}=n_1\stackrel{}{c}_1+n_2\stackrel{}{c}_2(n_1,n_2=0,\pm 1,\pm 2,..)`$ where $`\stackrel{}{c}_1=(1,0),\stackrel{}{c}_2=(\frac{1}{2},\frac{\sqrt{3}}{2})`$ are the basic vectors of the triangular lattice, $`J_\stackrel{}{a}=J\left(1\eta \stackrel{}{a}(\stackrel{}{u}_{\stackrel{}{n}+\stackrel{}{a}}\stackrel{}{u}_\stackrel{}{n})\right)=J\left(1\eta _{i,j}a_iu_{ij}a_j\right)`$ is the exchange integral ($`J>0`$ is the exchange constant), $`\stackrel{}{u}_\stackrel{}{n}`$ are the atom displacements from their equilibrium (without magnetic interaction), and $`\stackrel{}{a}`$ is the vector which connects a site with nearest neighbors, $`\eta =d\mathrm{ln}J/d|\stackrel{}{a}|`$. We consider a small spatially smooth deformation. Therefore it was convenient to introduce the components of the strain tensor $`u_{ij}`$: $`(\stackrel{}{u}_{\stackrel{}{n}+\stackrel{}{a}}\stackrel{}{u}_\stackrel{}{n})_i=\frac{1}{2}_j\left(\frac{u_i}{n_j}+\frac{u_j}{n_i}\right)a_j_ju_{ij}a_j,`$ ($`i,j=x,y).`$

The elastic energy of the two-dimensional triangular lattice in the continuum approximation coincides with the elastic energy of the isotropic plane:$`\mathrm{\Phi }=\frac{1}{2}K(u_{xx}+u_{yy})^2+\frac{1}{2}\mu \left((u_{xx}u_{yy})^2+4u_{xy}^2\right)p_xu_{xx}p_yu_{yy}`$ where $`K`$ and $`\mu `$ are the elastic modulus of the system and $`p_i(i=x,y)`$ are the component of the uniaxial pressure.

Let us use the transformation (see e.g. ) to the local frame of reference $`S_\stackrel{}{n}^x=\stackrel{~}{S}_\stackrel{}{n}^x\mathrm{cos}\stackrel{}{q}\stackrel{}{n}+\stackrel{~}{S}_\stackrel{}{n}^z\mathrm{sin}\stackrel{}{q}\stackrel{}{n},S_\stackrel{}{n}^z=\stackrel{~}{S}_\stackrel{}{n}^z\mathrm{cos}\stackrel{}{q}\stackrel{}{n}\stackrel{~}{S}_\stackrel{}{n}^x\mathrm{sin}\stackrel{}{q}\stackrel{}{n},S_\stackrel{}{n}^y=\stackrel{~}{S}_\stackrel{}{n}^y`$ in which quantization axis for the spins at each site coincide with its classical direction $`(\mathrm{cos}\stackrel{}{q}\stackrel{}{n},0,\mathrm{sin}\stackrel{}{q}\stackrel{}{n})`$ which is determined by the vector $`\stackrel{}{q}`$. By using the Holstein-Primakoff spin-representation $`\stackrel{~}{S}_\stackrel{}{n}^x=\sqrt{S/2}(b_\stackrel{}{n}^{}+b_\stackrel{}{n}),\stackrel{~}{S}_\stackrel{}{n}^y=i\sqrt{S/2}(b_\stackrel{}{n}^{}b_\stackrel{}{n}),\stackrel{~}{S}_\stackrel{}{n}^z=Sb_\stackrel{}{n}^{}b_\stackrel{}{n}`$ with $`b_\stackrel{}{n}^{}`$, $`(b_\stackrel{}{n})`$ being the creation (destruction) Bose-operator of the spin-excitation on the $`\stackrel{}{n}`$-th site and applying the Bogolyubov transformation $`b_\stackrel{}{n}=N^{1/2}e^{i\stackrel{}{k}\stackrel{}{n}}\left(\mathrm{cosh}(\theta _\stackrel{}{k})\alpha _\stackrel{}{k}+\mathrm{sinh}(\theta _\stackrel{}{k})\alpha _\stackrel{}{k}^{}\right)`$ for these operators ($`\stackrel{}{k}`$ is the wave vector belonging to the first Brillouin zone and $`N`$ is the number of atoms in the crystal), the Hamiltonian (1) in the spin-wave approximation can be represented in the form $`H=NE_{gr}+_\stackrel{}{k}\omega _\stackrel{}{k}\alpha _\stackrel{}{k}^{}\alpha _{veck}.`$ where

$`E_{gr}=E_{cl}{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k}}{}}\mathrm{sinh}^2(\theta _\stackrel{}{k})\omega _\stackrel{}{k}`$ (2)
is the ground state energy per atom in the spin-wave approximation. The coefficients of the Bogolyubov transformation are given by $`\mathrm{tanh}^2(\theta _\stackrel{}{k})=\left(f_\stackrel{}{k}g_\stackrel{}{k}\right)^2\left(f_\stackrel{}{k}+g_\stackrel{}{k}\right)^2`$ where

$`f_\stackrel{}{k}^2=S\left(2J(\stackrel{}{q})+J(\stackrel{}{k}\stackrel{}{q})+J(\stackrel{}{k}+\stackrel{}{q})\right),`$ (3)
$`g_\stackrel{}{k}^2=2S\left(J(\stackrel{}{k})J(\stackrel{}{q})\right)`$ (4)
with $`J(\stackrel{}{k})=\frac{1}{2}_\stackrel{}{a}J_\stackrel{}{a}\mathrm{cos}\stackrel{}{k}\stackrel{}{a}`$, and

$`\omega _\stackrel{}{k}=f_\stackrel{}{k}g_\stackrel{}{k}`$ (5)
is the magnon energy. $`E_{cl}=S^2J(\stackrel{}{q})`$ is the energy of the spin subsystem per atom in the classical approximation.

The on-site magnetization $`M\stackrel{~}{S}_\stackrel{}{n}^z=Sb_\stackrel{}{n}^{}b_\stackrel{}{n}`$ can be also obtained by using the Bogolyubov transformation. As a result we get

$`M=S{\displaystyle \frac{1}{N}}{\displaystyle \underset{\stackrel{}{k}}{}}\mathrm{sinh}^2(\theta _\stackrel{}{k}).`$ (6)
Our goal now is to minimize the total ground state energy of the system

$`=\mathrm{\Phi }+E_{gr}E_{gr}+J\left({\displaystyle \frac{\kappa }{2}}(ϵ_1^2+ϵ_2^2)pϵ_1\right)`$ (7)
with respect to the variational parameters which are the components of the deformation tensor $`ϵ_1=\eta (u_{xx}u_{yy})/2,ϵ_2=\eta u_{xy}`$ and the vector $`\stackrel{}{q}`$. Here the magneto-elastic constant $`\kappa =4\mu /\left(J\eta ^2\right)`$ characterizes the stiffness of the lattice in terms of the intensity of the spin-lattice interaction, $`p=(p_xp_y)/\left(\eta J\right)`$is the dimensionless uniaxial pressure. In Eq. (7) we omitted the terms of spin-lattice interaction which correspond to the dilatation (contraction) of the lattice ($`(u_{xx}+u_{yy})/2`$) without changing its symmetry. These terms don’t change the qualitative behavior of the system. The equations for $`\stackrel{}{q}`$ and $`ϵ_i`$ have the form

$`{\displaystyle \frac{}{\stackrel{}{q}}}=0,{\displaystyle \frac{}{ϵ_i}}=0`$ (8)
Let us consider first the case zero pressure: $`p=0`$. In the classical approximation when $`S\mathrm{}`$ the energy of the system can be represented in the form

$`{\displaystyle \frac{_{cl}}{J}}=S^2J(\stackrel{}{q})+{\displaystyle \frac{\kappa }{2}}(ϵ_1^2+ϵ_2^2)`$ (9)
¿From Eq.(9) we obtain that the set of equations (8) has three types of solutions:

i) H-phase: $`\stackrel{}{q}_H=4\pi /3(\mathrm{cos}\varphi ,\mathrm{sin}\varphi ),\varphi =0,\pm \pi /3.`$ The spin structure corresponding to each $`\stackrel{}{q}_H`$ is a three-sublattice antiferromaget.Three different $`\stackrel{}{q}_H`$ represent three possible antiferromagnetic domains. The lattice structure is an undistorted triangular lattice $`(ϵ_1=ϵ_2=0)`$ with the group symmetry $`D_{6h}`$.

ii) L-phase: $`\stackrel{}{q}_L=2\pi /\sqrt{3}(\mathrm{sin}\varphi ,\mathrm{cos}\varphi ),\varphi =0,\pm \pi /3.`$ The spin structure corresponding to each $`\stackrel{}{q}_L`$ is a two-sublattice antiferromagnet. The lattice structure for $`\varphi =0`$ is determined by the deformation tensors $`ϵ_1=2S^2/\kappa ,ϵ_2=0`$ while for two other $`\stackrel{}{q}_L`$ the deformation tensors can be obtained by rotations by $`\pm \pi /3`$. The lattice is characterized by the group symmetry $`D_{2h}`$.

iii) S-phase: $`\stackrel{}{q}_s=(q_s,2\pi /\sqrt{3}),ϵ_1=\kappa ^1\left(\mathrm{cos}q_s+\mathrm{cos}(q_s/2)\right),ϵ_2=0`$ where $`\mathrm{cos}(q_s/2)=\left(\sqrt{9+32\kappa /S^2}5\right)/8`$. Two more vectors $`\stackrel{}{q}`$ and deformation tensors $`ϵ_1`$ and $`ϵ_2`$ can be obtained from the above equation by rotations by $`\pm \pi /3`$. The spin structure is a spiral antiferromagnet.

We consider the stability of obtained solutions for the case when

$`ϵ_10,ϵ_2=0,0q_x{\displaystyle \frac{2\pi }{3}},q_y={\displaystyle \frac{2\pi }{\sqrt{3}}},`$ (10)
In the subspace (10) the stability condition reads

$`{\displaystyle \frac{^2}{ϵ_1^2}}>0,{\displaystyle \frac{^2}{ϵ_1^2}}{\displaystyle \frac{^2}{q_x^2}}\left({\displaystyle \frac{^2}{q_xϵ_1}}\right)^2>\mathrm{\hspace{0.17em}0}`$ (11)
One can obtain from Eqs (9) and (11) that

* the $`L`$-phase is stable when $`\kappa /S^2<\frac{9}{4}`$
* both $`H`$\- and $`L`$-phases are stable when $`\frac{9}{4}<\kappa /S^2<5.`$ For $`\kappa /S^2=4`$ the H- and L-phases have the same energy while for $`\kappa /S^2<4`$ the L-phase becomes energetically more favorable.
* the $`H`$-phase is stable when $`\kappa /S^2>5`$.

The spiral S-phase does not correspond to a minimum of the ground state energy. In the subspace (10) the spiral phase for $`\frac{9}{4}<\kappa /S^2<5`$ corresponds to the saddle point which separates two stable phases $`H`$ and $`L`$.

The wave vectors $`\stackrel{}{q}_H`$ and $`\stackrel{}{q}_L`$ determine also the magnetic structure of the H- and L-phases in the spin-wave approximation. This is not true, however, for the spiral phase where quantum fluctuations change the value of the spiral wave-vector $`q_s`$.

Evaluating numerically the integrals in the r.h.s. of Eq. (2) for $`\stackrel{}{q}=\stackrel{}{q}_H,ϵ_1=0`$ one can obtain that the ground state energy of the symmetric triangular phase in the spin-wave approximation is $`_HJ\left(\frac{3}{2}S(S+1)+1.1723S\right)`$. We have also found that the H-phase is stable for $`\kappa >\kappa _H`$ where

$$\frac{\kappa _H}{S^2}=\{\begin{array}{cc}\hfill 6.86,& \text{for }S=\frac{1}{2},\hfill \\ \hfill 3.85,& \text{for }S=1,\hfill \\ \hfill 3.26,& \text{for }S=\frac{3}{2},\hfill \\ \hfill 2.99,& \text{for }S=2\hfill \end{array}$$
(12)
Considering the low-symmetry L-phase which corresponds to the wave vector $`\stackrel{}{q}_L=2\pi (0,1/\sqrt{3})`$ we obtain from Eqs (3), (5) that its energy in the spin-wave approximation is determined by the expression

$`{\displaystyle \frac{_L}{J}}=S(S+1)(1+2ϵ_1)+{\displaystyle \frac{\kappa ϵ_1^2}{2}}+SI(ϵ_1)`$ (13)
$`I(ϵ_1)=(2+ϵ_1){\displaystyle \frac{4}{\pi ^2}}{\displaystyle \underset{0}{\overset{\frac{\pi }{2}}{}}}𝑑y(12{\displaystyle \frac{1ϵ_1}{2+ϵ_1}}\mathrm{sin}^2y)E(m)`$ (14)
where $`E(m)`$ is the elliptic integral of the second kind  with the parameter $`m=\mathrm{cos}^2y\left(12\frac{1ϵ_1}{2+ϵ_1}\mathrm{sin}^2y\right)^2`$. One can obtain from Eq. (5) that in the $`\stackrel{}{k}0`$ limit the magnon dispersion in the L-phase has the form $`\omega _L(\stackrel{}{k})=2JS(2+ϵ_1)\sqrt{(5ϵ_12)k_x^2+3(2+ϵ_1)k_y^2}`$. This means that the magnons in the low-symmetry L-phase can exist only for $`ϵ_1>0.4`$ and only in this interval we may look for extrema of the function $`_L`$. It is seen from Eqs (1) that the quantity $`J_1=J(1+ϵ_1/2)`$ $`\left(J_2=J(1ϵ_1)\right)`$ corresponds to the antiferromagnetic exchange between first (second) neighbors in the spin-wave theory of Heisenberg antiferromagnet on an anisotropic triangular lattice . Therefore the interval $`1>ϵ_1>0.4`$ corresponds to the interval $`0<J_2/J_1<\mathrm{\hspace{0.17em}1}/2`$ in the spin-wave theory of Heisenberg antiferromagnet on an anisotropic triangular lattice.

In Fig. 1 we plot the deformation tensor $`ϵ_1`$ which provides an extremum of $`_L`$, as a function of the coupling parameter $`\kappa `$ for three spin values $`S=1/2,1,3/2`$. It is seen that the equation $`d/dϵ_1=0`$ has solution $`ϵ_1`$only for $`\kappa <\kappa _L`$. The critical value $`\kappa _L`$ which determines the interval of the existence of the L-phase depends on the spin $`S`$. It increases when the spin $`S`$ increases, e.g.

$$\frac{\kappa _L}{S^2}=\{\begin{array}{cc}\hfill 3.97,& \text{for }S=\frac{1}{2},\hfill \\ \hfill 4.47,& \text{for }S=1,\hfill \\ \hfill 4.64,& \text{for }S=\frac{3}{2},\hfill \\ \hfill 4.73,& \text{for }S=2\hfill \end{array}$$
(15)
but even for $`S=4`$ it is still less ($`\kappa _L/S^24.75`$) than the critical value for this parameter obtained in the classical approximation ($`\kappa _L/S^2=5`$). The reason for this is the strong contribution to the effective elastic energy of the system from quantum fluctuations (the second term in the r.h.s. of Eq. (13)).The quantum fluctuations change both the stiffness of the lattice (the second derivative of the free energy (13) with respect to the deformation $`ϵ_1`$) and the constant of the spin-elastic interaction ( the first derivative of the free energy). A single-valued monotonic dependence $`ϵ_1(\kappa )`$ is obtained for $`S1`$. For $`S=1/2`$ the dependence becomes multi-valued, i.e. there is an interval of $`\kappa `$ where two values $`ϵ_1^{(1)}`$ and $`ϵ_1^{(2)}`$ ($`ϵ_1^{(1)}<ϵ_1^{(2)}`$ ) of the deformation tensor correspond to each value of the coupling constant $`\kappa `$. Only the solution which corresponds to the deformation $`ϵ_1^{(2)}`$ corresponds to the minimum of the effective elastic energy (13).

To find the stability region of the collinear L-phase we numerically evaluated integrals in the l.h.s. of the inequality (11) and found that it holds in the interval $`\kappa <\kappa _L`$. Comparing the stability limits for H- and L-phases give in Eqs (12) and (15) we find the following interesting results. For $`S\mathrm{\hspace{0.17em}1}`$ we have always $`\kappa _H<\kappa _L`$ like in the classical approach where the stability regions for H- and L-phases overlap. However, for $`S=1/2`$ this is not the case since there is an interval $`\kappa _L<\kappa <\kappa _H`$ where neither the L- nor the H-phase exist.

The nature of the phase in the interval $`(\kappa _L,\kappa _H)`$ or in other words for $`ϵ_1<0.4`$, cannot be clarified in the framework of the spin-wave approach. The reason why this approach fails is the following. As it was mentioned above the vectors $`\stackrel{}{q}_s`$ obtained in the classical approach are not solutions of the extrema conditions (8) in the spin-wave approximation. On the other hand, as it is seen from Eqs (5) and (3) the magnon frequency $`\omega _\stackrel{}{k}`$ is real in the interval $`ϵ_1<0.4`$ only for $`\stackrel{}{q}`$ obtained in the classical approach. For other $`\stackrel{}{q}`$ values $`\omega _\stackrel{}{k}^2`$ is not positive definite. Thus to check the stability of the corresponding phase is stable one should calculate the magnon dispersion beyond the spin-wave approach.

Let us consider how the L-phase evolves in the presence of an uniaxial pressure $`p`$. From Eqs (7) and (8) we get

$`\stackrel{}{q}={\displaystyle \frac{2\pi }{\sqrt{3}}}(0,1),ϵ_2=0,`$ (16)
$`p=\kappa ϵ_12S(S+1)+S{\displaystyle \frac{dI(ϵ_1)}{dϵ_1}}`$ (17)
In Fig. 2 we show results for the on-site magnetization in the L-phase for $`S=\frac{1}{2}`$ as a function of the uniaxial pressure $`p`$. Fig. 2 and in particular its inset shows that our results are qualitatively different from those of a recent spatially anisotropic triangular lattice study . In contrast to the $`(J_1,J_2)`$ -model  where the ratio $`J_2/J_1`$ was a free parameter in our model the corresponding parameter in a compressible triangular lattice $`J_2/J_1=2(1ϵ_1)/(2+ϵ_1)`$ is determined self-consistently from the extrema conditions (8). It is seen that in the interval of the existence of the L-phase we have $`J_2/J_1<1/2`$ $`(ϵ_1>0.4`$) and the on-site magnetization is finite. It is interesting that there exists an interval where two values of the magnetization $`M`$ correspond to each value of the pressure $`p`$. One can show, however, that only the state with a larger value of $`M`$ is stable.

We conclude that the Heisenberg antiferromagnet on an compressible triangular lattice differs qualitatively from the Heisenberg antiferromagnet on an anisotropic triangular lattice. Self-consistent treatment of the elastic degrees of freedom shows that the spiral magnetic phase is unstable in the classical approximation. This approach shows also that in the spin-wave approximation the interaction between quantum fluctuations and elastic degrees of freedom stabilizes the collinear L-phase.

Acknowledgements

Yu. Gaididei is grateful for the hospitality of the University of Bayreuth

where this work was performed. Partial support from the DRL grant Nr.: UKR-002-99 is also acknowledged.


text:
# 1 The sensitivity range reached by the planned GENIUS experiment in the mass parameter 𝑚 and the mixing angle variables 𝑏=sin²{𝜃₁₂}+sin²{𝜃₁₃} and 𝑏'=sin²{𝜃₁₂}-sin²{𝜃₁₃} in Model A. The upper and/or lower bounds of ±0.016 for 𝑏 and 𝑏' are obtained from Bugey experiment []. The area bounded by the dashed line illustrates Cases I (𝑏 as vertical axis) and III (𝑏' as vertical axis) (see eq.()) whereas the area bounded by solid line illustrates the sensitivity range in Cases II (𝑏) and IV (𝑏'). Cases V-VIII are related to Cases I-IV in a way described in the text.

HIP-2000-11/TH

Neutrinoless double beta decay

in four-neutrino models

Anna Kalliomäki 111E-mail address: amkallio@pcu.helsinki.fi and Jukka Maalampi 222E-mail address: maalampi@pcu.helsinki.fi

Department of Physics, Theoretical Physics Division

University of Helsinki, Finland

ABSTRACT

The most stringent constraint on the so-called effective electron neutrino mass from the present neutrinoless double beta decay experiments is $`|M_{ee}|<0.2`$ eV, while the planned next generation experiment GENIUS is anticipated to reach a considerably more stringent limit $`|M_{ee}|<0.001`$ eV. We investigate the constraints these bounds set on the neutrino masses and mixings of neutrinos in four-neutrino models where there exists a sterile neutrino along with the three ordinary neutrinos. We find that the GENIUS experiment would be sensitive to the electron neutrino masses down to the limit $`m_{\nu _e}<\mathrm{\hspace{0.33em}0.024}`$ eV in such a scenario.

There exist strong experimental indications that neutrinos have mass and that they mix. The results of the Super-Kamiokande experiment  on neutrinos produced in Earth’s atmosphere by cosmic rays provide a convincing evidence of neutrino oscillations, and the results of the solar neutrino experiments  point to the same direction. Results from the laboratory experiment by the LSND collaboration  may also be an indication of neutrino oscillations but they still await a future confirmation, e.g. at the MiniBoone and I216 experiments , in particular so because measurements in KARMEN detector  exclude most of the region in the parameter space favoured by the LSND.

Neutrino oscillation probabilities are determined by squared mass differences $`\delta m_{ji}^2=m_j^2m_i^2`$, where $`m_i`$ is the mass of the massive neutrino state $`\nu _i`$, and the elements of the neutrino mixing matrix $`U`$ that connects the massive neutrino states $`\nu _i`$ and the flavour neutrino states $`\nu _\alpha `$ through the relation $`\nu _\alpha =_iU_{\alpha i}\nu _i`$. For explaining the solar and atmospheric neutrino data two mass-squared difference scales, $`\delta m_{\mathrm{atm}}^210^310^2\mathrm{eV}^2`$ for atmospheric neutrinos and $`\delta m_{\mathrm{sun}}^210^5\mathrm{eV}^2`$ or $`10^{10}\mathrm{eV}^2`$ for solar neutrinos, are needed, and the data can be described within a three-neutrino model. If the LSND result is not disregarded, a three-neutrino framework is not sufficient but one has to go beyond it and consider a model with four (or more) neutrinos as three different levels of mass-squared differences are needed, the third one being $`\delta m_{\mathrm{LSND}}^2=0.32`$ eV2. In addition to the three known neutrino flavours $`\nu _e,\nu _\mu `$ and $`\nu _\tau `$ one has to incorporate at least one light neutrino $`\nu _s`$ that is sterile, i.e. it lacks Standard Model (SM) gauge interactions .

As was pointed out in ref. , the only four-neutrino mass patterns that are consistent with all oscillation data are those where there are two pairs of neutrinos with close masses, the ”sun-pair” (mass states $`\nu _0`$ and $`\nu _1`$) and ”atm-pair” (mass states $`\nu _2`$ and $`\nu _3`$) corresponding to the squared-mass differences $`\delta m_{\mathrm{atm}}^2`$ and $`\delta m_{\mathrm{sun}}^2`$, respectively, separated by a gap corresponding to $`\delta m_{\mathrm{LSND}}^2`$ .

Another oscillation constraint follows from cosmology. Depending on the value of $`\delta m^2`$, a large active-neutrino mixing could bring sterile neutrinos into thermal equilibrium thereby increasing the effective number of light neutrinos , $`N_\nu `$, which would affect the big bang nucleosynthesis (BBN). This leads to quite tight constraints on the active-sterile mixings .

In the case neutrinos are Majorana particles , which is quite likely from the theoretical point of view, another useful source of experimental information on neutrino masses and mixings is provided by the neutrinoless double beta ($`0\nu \beta \beta `$) decay. The $`0\nu \beta \beta `$ decay is sensitive to the values of the neutrino masses themselves, in contrast with the oscillation data which depend only on the squared-mass differences. Also, it provides information on the Majorana CP-phases of neutrinos, to which the oscillation phenomena are practically insensitive. In the present work we shall study the constraints the $`0\nu \beta \beta `$ decay places on the neutrino masses $`m_i`$, the neutrino mixing matrix $`U`$ and the relative CP-phases of neutrinos in four-neutrino models consisting of three active neutrinos $`\nu _e,\nu _\mu `$ and $`\nu _\tau `$ and a sterile neutrino $`\nu _s`$ with the 2+2 mass pattern described above.

The quantity probed in the $`0\nu \beta \beta `$ decay experiments is the absolute value of the effective mass

$$M_{ee}=\underset{i}{}m_iU_{ei}^2.$$
(1)
The most recent experimental upper bound for it from the Moscow-Heidelberg experiment is (at 90% C.L.)

$$|M_{ee}|<0.2\mathrm{eV}.$$
(2)
In turns out that this bound is in the case of four-neutrino models less restrictive than in the case of the tree-neutrino models . A substantial strengthening of the present bound is, however, foreseen in the future. In the planned GENIUS experiment one expects to reach after the first year of operation the limit 0.01 eV and eventually the limit

$$|M_{ee}|<0.001\mathrm{eV}(\mathrm{GENIUS}).$$
(3)
We will study the restrictions this bound would set on the parameters of four-neutrino models, in particular on the mass of the electron neutrino $`\nu _e`$ and the mixing of $`\nu _e`$ with $`\nu _\mu `$ and $`\nu _\tau `$.

The constraint from the $`0\nu \beta \beta `$ decay is sensitive to the relative CP-phases $`\eta _i`$ of the mass states $`\nu _i`$. We will concentrate in the CP-conserving case where $`\eta _i=\pm 1`$. In this case one can write (1) in the form

$$M_{ee}=\underset{i}{}\eta _im_i|U_{ei}|^2.$$
(4)
The most general mixing matrix $`U`$ that describes the connection between the flavour states $`\nu _s,\nu _e,\nu _\mu ,\nu _\tau `$ and the Majorana mass eigenstates $`\nu _0,\nu _1,\nu _2,\nu _3`$ can be parametrized in terms of $`6`$ rotation angles and $`6`$ phases as follows (see, e.g., ref. ):

$$U=\left(\begin{array}{cccc}c_{01}c_{02}c_{03}& c_{02}c_{03}s_{01}^{}& c_{03}s_{02}^{}& s_{03}^{}\\ & & & \\ c_{01}c_{02}s_{03}s_{13}^{}& c_{02}s_{01}^{}s_{03}s_{13}^{}& s_{02}^{}s_{03}s_{13}^{}& c_{03}s_{13}^{}\\ c_{01}c_{13}s_{02}s_{12}^{}& c_{13}s_{01}^{}s_{02}s_{12}^{}& +c_{02}c_{13}s_{12}^{}& \\ c_{12}c_{13}s_{01}& +c_{01}c_{12}c_{13}& & \\ & & & \\ c_{01}c_{02}c_{13}s_{03}s_{23}^{}& c_{02}c_{13}s_{01}^{}s_{03}s_{23}^{}& c_{13}s_{02}^{}s_{03}s_{23}^{}& c_{03}c_{13}s_{23}^{}\\ +c_{01}s_{02}s_{12}^{}s_{13}s_{23}^{}& +s_{01}^{}s_{02}s_{12}^{}s_{13}s_{23}^{}& c_{02}s_{12}^{}s_{13}s_{23}^{}& \\ c_{01}c_{12}c_{23}s_{02}& c_{12}c_{23}s_{01}^{}s_{02}& +c_{02}c_{12}c_{23}& \\ +c_{12}s_{01}s_{13}s_{23}^{}& c_{01}c_{12}s_{13}s_{23}^{}& & \\ +c_{23}s_{01}s_{12}& c_{01}c_{23}s_{12}& & \\ & & & \\ c_{01}c_{02}c_{13}c_{23}s_{03}& c_{02}c_{13}c_{23}s_{01}^{}s_{03}& c_{13}c_{23}s_{02}^{}s_{03}& c_{03}c_{13}c_{23}\\ +c_{01}c_{23}s_{02}s_{12}^{}s_{13}& +c_{23}s_{01}^{}s_{02}s_{12}^{}s_{13}& c_{02}c_{23}s_{12}^{}s_{13}& \\ +c_{01}c_{12}s_{02}s_{23}& +c_{12}s_{01}^{}s_{02}s_{23}& c_{02}c_{12}s_{23}& \\ +c_{12}c_{23}s_{01}s_{13}& c_{01}c_{12}c_{23}s_{13}& & \\ s_{01}s_{12}s_{23}& +c_{01}s_{12}s_{23}& & \end{array}\right),$$
(5)
where $`c_{jk}\mathrm{cos}\theta _{jk}`$ and $`s_{jk}\mathrm{sin}\theta _{jk}e^{i\delta _{jk}}`$. For our discussion relevant is the second row which gives the composition of the electron neutrino in terms of the mass eigenstates neutrinos. It depends on several mixing parameters but fortunately simplifies considerably when the existing constraints from laboratory experiments and cosmology are taken into account.

Assuming the 2+2 mass hierarchy among neutrino masses, indicating $`\delta m_{02}^2\delta m_{03}^2\delta m_{12}^2\delta m_{13}^2\delta m_{LSND}^2`$, the Bugey short baseline experiment , measuring the probability $`P(\nu _e\nu _e)1A^{ee}\mathrm{sin}^2(L\delta m_{\mathrm{LSND}}^2/4E)`$, yields the bound

$$A^{ee}=4(|U_{e2}|^2+|U_{e3}|^2)(1|U_{e2}|^2|U_{e3}|^2)<0.06,$$
(6)
or

$$|c_{02}c_{13}s_{12}^{}s_{02}^{}s_{03}s_{13}^{}|^2+|c_{03}s_{13}^{}|^2<0.016.$$
(7)
The angles $`\theta _{02}`$ and $`\theta _{03}`$ are further constrained by the big bang nucleosynthesis (BBN) as large mixings would, as mentioned above, increase the number of the effective degrees of freedom $`N_\nu `$ by bringing the sterile neutrino into thermal equilibrium. According to a recent analysis , the two-flavour $`\nu _{\tau ,\mu }\nu _s`$ mixings should obey the constraint (valid for $`|\delta m^2|<2.510^3\mathrm{eV}^2`$)

$$|\delta m^2|\mathrm{sin}^2\theta <710^5\mathrm{eV}^2$$
(8)
in order not to lead to a conflict with the BBN limit $`N_\nu <\mathrm{\hspace{0.33em}3.2}`$ . Since now $`\delta m^2=\delta m_{\mathrm{LSND}}^2=0.32`$ eV2, the mixing angle $`\theta `$ should be $`<\mathrm{\hspace{0.33em}10}^2`$, and it is therefore conceivable to assume $`s_{02},s_{03}1`$. With this approximation one derives from (7) the bounds

$`b`$ $``$ $`\mathrm{sin}^2\theta _{12}+\mathrm{sin}^2\theta _{13}<\mathrm{\hspace{0.33em}0.016},`$
$`|b^{}|`$ $``$ $`|\mathrm{sin}^2\theta _{12}\mathrm{sin}^2\theta _{13}|<\mathrm{\hspace{0.33em}0.016},`$ (9)
and the following approximative composition of the electron neutrino in terms of the mass eigenstates:

$$\nu _e=\mathrm{sin}\theta _{01}e^{i\delta _{01}}\nu _0\mathrm{cos}\theta _{01}\nu _1+\mathrm{sin}\theta _{12}e^{i\delta _{12}}\nu _2+\mathrm{sin}\theta _{13}e^{i\delta _{13}}\nu _3.$$
(10)
Let us note that the cosmological argument used above remains still somewhat controversial. Some analysis allow for the number of the relativistic degrees of freedom, instead of the bound $`N_\nu <\mathrm{\hspace{0.33em}3.2}`$ quoted above, values close to four , in which case BBN would not place any constraint on the angles $`s_{02},s_{03}`$. The difference is due to conflicting results on the primordial deuterium abundance. Furthermore, the BBN constraint depends on the lepton asymmetry in the early universe, which might considerably weaken the bounds  (in the case the predominantly sterile mass eigenstate is lighter than its active mixing partner). Anyhow, we will assume in the following that the sterile neutrino mixes considerably only with the electron neutrino and consequently that the approximations (9) and (10) are valid.

In this approximation the solar neutrino oscillation is described entirely in terms of the angle $`\theta _{01}`$ associated with the $`\nu _e\nu _s`$ mixing. According to an overall analysis of all solar neutrino data, only the small angle matter solution (SAMSW) with $`210^3<\mathrm{sin}^22\theta _{01}<\mathrm{\hspace{0.33em}10}^2`$ leads to a reasonable fit in the active-sterile mixing case .

There are two mass patterns consistent with the oscillation data, differing in as to which one of the pairs, the sun-pair or the atm-pair, is lighter. We will call Model A the pattern where the sun-pair is the light one and the atm-pair is the heavy one, and the case where the mass order is the opposite will be called Model B. The mass scale of the lighter neutrino pair is denoted by $`m`$, and the mass gap between the pairs is denoted by $`\mathrm{\Delta }m`$, so that the mass scale of the heavier pair is $`m+\mathrm{\Delta }m`$. We will neglect the mass differencies inside the two pairs by assuming they to be small as compared with the mass gap $`\mathrm{\Delta }m`$.

Studies of the anisotropies of the cosmic microwave background radiation and the large scale structure of the universe provide information on the sum of neutrino masses, $`_im_{\nu _i}4m+2\mathrm{\Delta }m`$. The contribution of neutrinos on the dark matter content of the Universe is in units of the critical density given by $`\mathrm{\Omega }_\nu =h^2_im_{\nu _i}/94`$ eV, where $`h`$ is the dimensionless Hubble constant ($`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$). In the previously popular hot+cold dark matter scenario  the neutrino contribution to the energy density of the universe is assumed to be $`\mathrm{\Omega }_\nu 0.2`$, but the more recent observations, indicating the existence of a large cosmological constant, imply that this should be taken as a conservative upper limit only . With the recent measured value of the Hubble constant , $`h0.71`$, this leads to the upper limit $`_im_{\nu _i}<\mathrm{\hspace{0.33em}10}`$ eV. Recently it has been claimed that this limit improves to

$$\underset{i}{}m_{\nu _i}<\mathrm{\hspace{0.33em}5.5}\mathrm{eV}$$
(11)
when the information from the Ly$`\alpha `$ forest in quasar spectra is taken into account .

The range of the squared mass difference, $`\delta m_{\mathrm{LSND}}^2=0.32`$ eV2, indicated by the LSND data and allowed by the Bugey and the other short baseline experiments yields the following upper and lower bounds for the mass difference $`\mathrm{\Delta }m`$ as a function of the mass $`m`$:

$$m+\sqrt{m^2+0.3\mathrm{eV}^2}<\mathrm{\Delta }m<m+\sqrt{m^2+2\mathrm{eV}^2}.$$
(12)
Combining this with the cosmological bound (11) leads to the following bounds on $`m`$ and $`\mathrm{\Delta }m`$:

$`m`$ $`<`$$``$ $`1.3\mathrm{eV},`$
$`\mathrm{\Delta }m`$ $``$ $`{\displaystyle \frac{\delta m_{\mathrm{LSND}}}{_im_i/2}}\mathrm{\hspace{0.33em}0.1}\mathrm{eV}.`$ (13)
In the extreme case of $`m=1.3`$ eV and $`\mathrm{\Delta }m=0.1`$ eV the mass spectrum is quite degenerate with $`m_0m_11.3`$ eV and $`m_2m_31.4`$ eV (Model A) or $`m_0m_11.4`$ eV and $`m_2m_31.3`$ eV (Model B). At the other extreme of $`m=0`$, i.e. when the lighter neutrino pair is massless, the heavier neutrino pair would have its mass in the range 0.5 eV to 1.4 eV.

Let us proceed to examine the consequences of the constraints (2) and (3) on the neutrino masses and mixing angles $`\theta _{01},\theta _{12}`$ and $`\theta _{13}`$ for different CP-phase patterns $`(\eta _0,\eta _1,\eta _2,\eta _3)`$ in the two models A and B.

Model A. In this model the neutrino pair ($`\nu _0,\nu _1`$) responsible for the solar anomaly is lighter than the pair ($`\nu _2,\nu _3`$) responsible for the atmospheric neutrino anomaly, i.e. the mass spectrum is of the form

$`m_0m_1=m`$
$`m_2m_3=m+\mathrm{\Delta }m.`$ (14)
The $`0\nu \beta \beta `$ condition takes in this case the form

$$|m(\eta _0\mathrm{sin}^2\theta _{01}+\eta _1\mathrm{cos}^2\theta _{01})+(m+\mathrm{\Delta }m)(\eta _2\mathrm{sin}^2\theta _{12}+\eta _3\mathrm{sin}^2\theta _{13})|<a,$$
(15)
where for the present experimental bound $`a=0.2`$ eV and for the anticipated future bound of the GENIUS experiment $`a=0.001`$ eV. In order to find out the consequences of this condition one should examine separately the following different patterns of the relative CP-phases $`\eta =(\eta _0,\eta _1,\eta _2,\eta _3)`$:

$`\mathrm{I},\mathrm{II}`$ $`\eta =(1,1,\pm 1,\pm 1),`$ (16)
$`\mathrm{III},\mathrm{IV}`$ $`\eta =(1,1,\pm 1,1),`$
$`\mathrm{V},\mathrm{VI}`$ $`\eta =(1,1,\pm 1,\pm 1),`$
$`\mathrm{VII},\mathrm{VIII}`$ $`\eta =(1,1,\pm 1,1).`$
In the first four cases (I-IV), where $`\eta _0=\eta _1`$, the $`0\nu \beta \beta `$ condition is not affected by the solar neutrino mixing angle $`\theta _{01}`$.

In Cases I and II the $`0\nu \beta \beta `$ condition (15) becomes

$$|m\pm (m+\mathrm{\Delta }m)b|<a$$
(17)
where the plus-sign corresponds to Case I and the minus-sign to Case II, and $`b=\mathrm{sin}^2\theta _{12}+\mathrm{sin}^2\theta _{13}`$, as defined in eq. (9).

Case I. The LSND and other short baseline experiments constrain the value of the mass gap $`\mathrm{\Delta }m`$ between the two neutrino pairs within the range given in eq. (12). Combining this with the cosmological lower bound (13) yields for $`\mathrm{\Delta }m`$ the allowed range from 0.1 eV to 1.4 eV. One then concludes from (17) that with the present experimental sensitivity the $`0\nu \beta \beta `$ result does not set any constraint on the mixing angles $`\theta _{12},\theta _{13}`$ as $`b\mathrm{\Delta }ma`$. The experimental constraint implies just an upper limit for the mass scale $`m`$, $`m<(ab\mathrm{\Delta }m)/(1+b)a=0.2`$ eV.

The situation would be different in the case of the anticipated GENIUS bound $`a=0.001`$ eV. The condition (17) to be satisfied one must now have $`ba/\sqrt{0.3}\mathrm{eV}2a/\mathrm{eV}=0.002`$ , which would mean tightening of the present bound (9) from the Bugey reactor experiment by one order of magnitude. Also the allowed values of $`m`$ are relatively small, $`m<a=0.001`$ eV.

Case II. The Case II does not differ from Case I in any essential way as far as the present $`0\nu \beta \beta `$ bound is concerned. With the GENIUS bound the upper limit of the mixing angles $`\theta _{12},\theta _{13}`$ is as stringest $`b<a/\sqrt{0.3}\mathrm{eV}2a/\mathrm{eV}=0.002`$ obtained for $`m=0`$. The situation is different compared with Case I in that now the condition (17) imposes also a non-zero lower limit on $`b`$ when $`m>a`$, otherwise there would not be sufficient cancellations among the varios terms to make the condition (15) to be valid. The quantity $`b`$ would be constrained into the range $`(ma)/\sqrt{m^2+2\mathrm{eV}^2}<b<(m+a)/\sqrt{m^2+0.3\mathrm{eV}^2}`$, where the upper limit is overtaken by the upper limit $`b<\mathrm{\hspace{0.33em}0.016}`$ from the short baseline experiments, when $`m>\mathrm{\hspace{0.33em}0.008}`$ eV. The sensitivity range of the GENIUS experiment in the $`(m,b)`$ parameter plane is presented in Fig. 1 (solid line). As can be read from the figure, the GENIUS experiment would in this case be sensitive to the values of $`m`$ down to 0.024 eV.

Cases III and IV. In Cases III and IV the $`0\nu \beta \beta `$ condition (15) becomes

$$|m\pm (m+\mathrm{\Delta }m)b^{}|<a,$$
(18)
where the plus-sign corresponds to Case III and the minus-sign to Case IV, and $`b^{}=\mathrm{sin}^2\theta _{12}\mathrm{sin}^2\theta _{13}`$, as defined in (9). Just like in Cases I and II, the present $`0\nu \beta \beta `$ bound does not probe the mixing angles $`\theta _{12}`$ and $`\theta _{13}`$ but sets a constraint merely on the mass $`m`$. Again, the GENIUS experiment would be sensitive to the values of these angles. In Case III the condition (18) restricts the allowed values of $`b^{}=\mathrm{sin}^2\theta _{12}\mathrm{sin}^2\theta _{13}`$, whose value is by other constraints limited to $`0.016<b^{}<\mathrm{\hspace{0.33em}0.016}`$, so that for example for $`m=0`$ only the values $`|b^{}|a/\sqrt{0.3}\mathrm{eV}2a/\text{eV}=0.002`$ are allowed. The condition (18) is never satisfied for the mass values $`m>\mathrm{\hspace{0.33em}0.024}`$ eV. This would be the sensitivity limit of the GENIUS on the electron neutrino mass also in this case. The allowed values of $`b^{}`$ are negative, $`(ma)/\sqrt{m^2+0.3\mathrm{eV}^2}<b^{}<(m+a)/\sqrt{m^2+2\mathrm{eV}^2}`$, except that for $`m<a`$ small positive values are also possible. The GENIUS sensitivity range for this case is presented in Fig. 1 with the dashed line. Case IV is a mirror image of Case III so that again for $`m=0`$ all values $`|b^{}|<a/\sqrt{0.3}`$ eV are possible, but now $`b^{}`$ can have mainly positive values: $`(ma)\sqrt{m^2+2\mathrm{eV}^2}<b^{}<(m+a)/\sqrt{m^2+0.3\mathrm{eV}^2}`$ (the solid line in Fig. 1). Note that this is almost the same region that was obtained in Case II for the quantity $`b`$. The limit that the GENIUS experiment could reach for the mass $`m`$ is again 0.024 eV.

Cases V and VI. In these cases the $`0\nu \beta \beta `$ condition depends also on the solar mixing angle $`\theta _{01}`$:

$$|m(\mathrm{sin}^2\theta _{01}\mathrm{cos}^2\theta _{01})\pm (m+\mathrm{\Delta }m)b|<a,$$
(19)
where the plus-sign corresponds to Case V and the minus-sign to Case VI. If only the small mixing angle (SMA) MSW solution is being considered, one can use the approximative value $`\mathrm{sin}^2\theta _{01}\mathrm{cos}^2\theta _{01}\pm 1`$, where the sign depends on whether $`\theta _{01}`$ is close to $`0`$ or $`\pi /2`$. If one assumes that angle $`\theta _{01}0,`$ condition (19) becomes

$$|m\pm (m+\mathrm{\Delta }m)b|<a.$$
(20)
As is apparent from eqs. (17) and (20), the $`0\nu \beta \beta `$ constraints in Cases II and V are in practical terms the same, and no separate treatment is needed for Case V. The same applies to Cases I and VI. On the other hand, if the angle $`\theta _{01}`$ is close to $`\pi /2`$, condition (20) must be replaced with

$$|m\pm (m+\mathrm{\Delta }m)b|<a.$$
(21)
and clearly Cases I and V as well as Cases II and VI correspond to each other as far as the $`0\nu \beta \beta `$ constraint is concerned.

It is noteworthy that the values of $`\theta _{01}`$ close to 0 and close to $`\pi /2`$ correspond here two physically distinct situations while in oscillation experiments, where transition probabilities depend on $`\theta _{01}`$ through $`\mathrm{sin}^22\theta _{01}`$, do not make any difference between the ranges $`0<\theta _{01}<\pi /4`$ and $`\pi /4<\theta _{01}<\pi /2`$.

Cases VII and VIII. These cases are practically identical to Cases III and IV. If $`\theta _{01}0,`$ the $`0\nu \beta \beta `$ condition (15) takes the form

$$|m\pm (m+\mathrm{\Delta }m)b^{}|<a,$$
(22)
where the plus-sign corresponds to Case VII and minus-sign to Case VIII. Comparison between eqs. (18) and (22) shows that in Cases IV and VII as well as in Cases III and VIII the $`0\nu \beta \beta `$ condition is practically the same. If $`\theta _{01}\pi /2,`$ the same applies to Cases III, VII and IV, VIII.

Model B. In this model the mass spectrum is of the form

$`m_0m_1=m+\mathrm{\Delta }m`$
$`m_2m_3=m.`$ (23)
Again we assume that mixing of the neutrino pair $`(\nu _0,\nu _1)`$ is mainly responsible for the solar neutrino anomaly and mixing of the pair $`(\nu _2,\nu _3)`$ for the atmospheric neutrino anomaly. The $`0\nu \beta \beta `$ condition becomes now

$$|(m+\mathrm{\Delta }m)(\eta _0\mathrm{sin}^2\theta _{01}+\eta _1\mathrm{cos}^2\theta _{01})+m(\eta _2\mathrm{sin}^2\theta _{12}+\eta _3\mathrm{sin}^2\theta _{13})|<a.$$
(24)
If LSND results are taken into account and one assumes that $`0.3\text{eV}^2<m_1^2m_3^2<\mathrm{\hspace{0.33em}2}\text{eV}^2`$, which yields $`m_1>\mathrm{\hspace{0.33em}0.5}\text{eV}`$, it is easy to see that eq. (24) is not realized for any value of $`m`$ . This is because the absolute values of $`b`$ and $`b^{}`$, defined in eq. (9), are much less than $`|\eta _0\mathrm{sin}^2\theta _{01}+\eta _1\mathrm{cos}^2\theta _{01}|1`$, so that eq. (24) can be written approximately as

$$|m+\mathrm{\Delta }m|<a,$$
(25)
which is untrue because $`|m+\mathrm{\Delta }m|=m_1`$ and even the present value of $`a`$ is less than the lower limit $`0.5`$eV of $`m_1`$. If we, on the other hand, disregard the LNSD results and let $`\mathrm{\Delta }m`$ be as low as $`0.1`$eV, with present value of $`a`$ we get some restrictions on the parameter values from eq. (24). Eight different $`\eta `$ combinations (16) are in principle possible also in this case but because (as was already mentioned) the term proportional to $`b`$ or $`b^{}`$ in eq. (24) is much smaller than the term proportional to angle $`\theta _{01}`$, it can safely be neglected. So, the approximative equation (25) is valid in every case, not depending on the choices of $`\eta _0`$ and $`\eta _1`$ or whether $`\theta _{01}`$ is $`0`$ or $`\pi /2`$. The allowed values for $`m`$ are in the range $`m<a\mathrm{\Delta }m`$.

We can summarize our results as follows. In the model A, i.e. when the neutrino pair (with mass $`m`$) responsible on the solar neutrino deficit is lighter than the neutrino pair (with mass $`m+\mathrm{\Delta }m`$) responsible on the atmospheric neutrino anomaly, the $`0\nu \beta \beta `$ constraint set by the future GENIUS experiment, $`M_{ee}<\mathrm{\hspace{0.33em}0.001}`$ eV, induces an upper limit for mass $`m`$ as a function of the quantity $`b=\mathrm{sin}^2\theta _{12}+\mathrm{sin}^2\theta _{13}`$ or the quantity $`b^{}=\mathrm{sin}^2\theta _{12}\mathrm{sin}^2\theta _{13}`$ depending on the pattern of the relative CP numbers of the neutrinos. There exists eight possible patterns of the relative CP numbers, Cases I to VIII, defined in eq. (16). In Case I, where all relative CP numbers are equal, and no cancellations between the contributions of different neutrinos hence occurs, the upper bound for $`m`$ is the stringest, $`m<\mathrm{\hspace{0.33em}0.001}`$ eV. For all the others Cases the absolute mass bound is about an order of magnitude less stringent, $`m<\mathrm{\hspace{0.33em}0.024}`$ eV. In Case I the value of $`b`$ must be small, $`b<a`$. In the other Cases the restriction on $`b`$ or $`b^{}`$ depend on the value of $`m`$ (see Fig. 1), and both a lower and an upper limit are obtained. For $`m<\mathrm{\hspace{0.33em}0.08}`$ eV the foreseen upper limit for $`b`$ or $`|b^{}|`$ is more stringent that the limit of the Bugey experiment. There is a crucial difference between the different relative CP number patterns as comes to the ensuing limitations on the mixing angles $`\theta _{12}`$ and $`\theta _{13}`$, as a bound on $`b`$ implies a bound on the both angles whereas a bound on $`b^{}`$ implies only a bound on the difference of the angles, not on each of the angles separately. Also, the smaller is the mass $`m`$, the smaller or the more degenerate are the angles $`\theta _{12}`$ and $`\theta _{13}`$, depending on whether the obtained bound is on $`b`$ or $`b^{}`$, respectively.

The Model B, i.e. when the neutrino pair responsible on the solar neutrino deficit is heavier than the neutrino pair responsible on the atmospheric neutrino anomaly, is in practical terms ruled out already by the present limit of $`M_{ee}<\mathrm{\hspace{0.33em}0.2}`$ eV.

This work has been supported by the Academy of Finland under the project no. 40677.