text:
# Nanodot-Cavity Electrodynamics and Photon Entanglement

## Abstract

Quantum electrodynamics of excitons in a cavity is shown to be relevant to quantum operations. We present a theory of an integrable solid-state quantum controlled-phase gate for generating entanglement of two photons using a coupled nanodot-microcavity-fiber structure. A conditional phase shift of $`O(\pi /10)`$ is calculated to be the consequence of the giant optical nonlinearity keyed by the excitons in the cavities. Structural design and active control, such as electromagnetic induced transparency and pulse shaping, optimize the quantum efficiency of the gate operation.

preprint: 031015-QCD

Semiconductor nanodot plays a key role in nanoscience as has been demonstrated by the electrical control of transport van der Wiel et al. (2003) and the optical control of quantum operations Li et al. (2003). Following the study of quantum electrodynamics of atoms in cavity (CQED) Berman (1994), effort is underway in the study of CQED of excitons in nanodots Wang et al. (2000). We report here the results of a theoretical study of excitons in CQED as illustrated by the proposal of a solid state controlled phase gate which entangles two photons.

Entangled photon pairs are the main stay of quantum information processing Bennett et al. (1996) and the controlled gate which conditions the dynamics of one photon on the state of the other also enables a key logic operation for quantum computation. There are two approaches to realize such gates: (1) linear optics with projective measurements Knill et al. (2001) and (2) nonlinear optics at the discrete photon level. The logic gate working with few-photon nonlinear optics requires impractical interaction length (e.g. several meters) in conventional Kerr media Sanders and Milburn (1989). To obtain giant optical nonlinearity for a two-photon logic gate, novel schemes have been demonstrated, e.g. the atom-cavity QED Turchette et al. (1995), or proposed, e.g. slow light in a coherently prepared atomic gas exhibiting electromagnetically induced transparency (EIT) Lukin and Imamo$`\overline{g}`$lu (2000).

The relevance of excitons in CQED is strengthened by the recent advances in solid state photonics and optoelectronics. We expect that the localization of the optical excitations would lead to ready integration of the solid state cavity devices with extant devices. Advances relevant to our proposal in semiconductor quantum devices include single photon sources operating at room temperature Michler et al. (2000); Lounis et al. (2000),high-Q microsheres and their coupling to nanodots Wang et al. (2000) and to fibers Spillane et al. (2003), and photonic lattice waveguides and cavities Scherer et al. (2002); Pang et al. (2003).

The qubit in our scheme is represented by two polarization states of a photon. In a quantum controlled phase gate, a two-photon state acquires a phaseshift conditional to their polarization configuration. The arrangement of our proposed device is given in Fig. 1(a). Two photons traveling along two optical fibers receive their interaction by coupling to two silicon microsphere cavities which are joined by a doped nanodot. The dot provides in theory Yao et al.  a strong third-order optical nonlinearity which is essential for a controlled interaction between two photons. Two cavities of different resonant frequencies are needed to afford control of coupling to either photon. They also act as an in-situ energy filter preventing two photons ending in the same fiber.

The photon scattering at the phase gate have inevitably some unwanted dynamics such as polarization-dependent reflection and motion-polarization entanglement. By relying on the transmission probability, the gate has the probabilistic nature as the linear optics procedure. The essential distinction lies in our use of the strong nonlinearity to provide definite interaction dynamics in the cavities versus the entanglement generated by the projective measurement. The probabilistic nature arises from the coupling of the photons to the solid state system, which is unavoidable in any system, but its effect can be ameliorated. Our solution is two pronged: to eliminate the linear reflection by EIT and to minimize motion-polarization entanglement by pulse shaping and system design.

The two $`\mathrm{LP}_{11}`$ modes in a step index optical fiber Cancellieri (1991) are chosen as the two polarization states $`|X`$ and $`|Y`$ for the qubits \[see Fig. 1(a)\]. The relevant modes in the microcavities are chosen to be the TE modes resonant with the nanodot transitions while the other TE modes and all TM modes are tuned far off-resonant for a small cavity ($`\mu `$m) Benner et al. (1980). The TE cavity mode can be excited only by an $`|X`$ photon in the fiber, whose coupling strength to the cavity on the left (right), $`\kappa _{L(R)}`$, are designed by adjusting the distance between the cavity and the fiber Spillane et al. (2003). Thus, only in the $`|XX`$ state do the two incoming photons interact via the cavity-dot coupling system, resulting a conditional phase-shift.

The strong photon-photon interaction induced by the dot-coupled cavities is favored by both the small cavity-mode volume and the large dipole moment of the nanodot transitions but the nanodot which contains a single active electron plays an essential role. The basic nonlinear optical process is illustrated with the aid of the energy structures in Figure 1(b). A strong magnetic field is applied along the $`x`$ direction to produce non-degenerate transitions from the electron spin states to the charged exciton states (trions), which are tuned respectively in resonance with the two cavity modes. The two split electron states are $`|x\pm (1/\sqrt{2})(e_+^{}\pm e_{}^{})|G`$, and the two degenerate trion states are $`|T\pm (1/\sqrt{2})(e_+^{}e_{}^{}h_{}^{}\pm e_+^{}e_{}^{}h_+^{})|G`$, where $`e_\pm ^{}`$ and $`h_\pm ^{}`$ create electron and hole spin states along the $`z`$ axis. The transition selection rules are: $`|x\pm |T`$ via the $`X`$ polarized field and $`|x\pm |T\pm `$ via the $`Y`$ polarized field. The spatial configuration of the cavity-dot structure is such that the TE modes are X-polarized at the site of the nanodot. The strong coupling between the trion state $`|T`$ and the cavity-dot state $`|x+,C_L`$ (or between $`|T+`$ and the cavity-dot state $`|x,C_R`$) mixes each pair into two split trion-polariton states, where $`C_{L(R)}`$ denotes the left (right) cavity mode. We denote the lower polariton states as $`|3`$ and $`|4`$, respectively. The four states, $`|x+`$, $`|x`$, $`|3`$, and $`|4`$, form the level structure for the optical nonlinearity and all other states are assumed far off resonance. This situation is well satisfied by the cavity-dot coupling $`g_{cav}0.5`$ meV, cavity linewidth $`0.1`$ meV, and Zeeman splitting $`g\mu _BB_x1`$ meV.

To induce an interaction between the photons from the left and right channels, a strong $`Y`$-polarized pump pulse is applied to resonantly couple the states $`|x`$ and $`|3`$. Consider the effect of this classical field in the dressed basis: $`|1d|x+|N`$, $`|3d|3|N`$, $`|2d|x|N+1`$, and $`|4d|4|N+1,`$ where the $`Y`$-polarized coherent field is approximated by the Fock state $`|N`$ with large $`N`$. Fig. 1(c) shows how the two $`X`$ photons on separate fibers which affect separately the modes in the left and right cavities are coupled by the $`Y`$ pump. The coupling strength $`\mathrm{\Omega }_c`$ between $`|3d`$ and $`|2d`$ is proportional to the electric field strength. Thus, the nonlinear optical coupling is readily manipulated by switching on and off the pump pulse. The classical pump pulse also increases the efficiency of the operation by cooling the spin system and by eliminating the linear reflection and absorption by laser cooling and EIT (see Eq. (1a)) Schmidt and Imamo$`\overline{g}`$lu (1996); Lukin and Imamo$`\overline{g}`$lu (2000).

The transformation of the polarization state of two photons is carried out by scattering theory. The initial state is specified by the density matrix $`\rho _i`$ in the basis set of the direct products of the polarization states, $`|\sigma _L\sigma _R`$ with $`\sigma =X`$ or $`Y`$, and of the wave vector states $`|k_L,k_R`$. The transmitted state $`\rho _f`$ is given by $`t\rho _it^{}`$, where $`t`$ is the transmission matrix. By design, $`t`$ is diagonal in the polarization states. The final density matrix of the two-photon polarizations is obtained by tracing $`\rho _f`$ over the wave vectors of the photons Mott and Massey (1965). The transmission $`t`$ is obtained via the $`T`$ matrix conserving the total energy. The linear and nonlinear scattering terms, illustrated by the Feynman diagrams, Fig. 1(d) and (e) indicating only the lowest order dressing terms by the $`Y`$ field, are non-perturbatively calculated Cohen-Tannoudji et al. (1992) as

$`T_{fi}^{\left(1\right)}=\delta _{\sigma _L,X}\delta _{k_R,k_R^{}}{\displaystyle \frac{\left|\kappa _L\right|^2}{2}}\times {\displaystyle \frac{E_{1d}+\mathrm{}ck_LE_{2d}}{\left(E_{1d}+\mathrm{}ck_LE_{3d}+i\mathrm{\Gamma }_3/2\right)\left(E_{1d}+\mathrm{}ck_LE_{2d}\right)\mathrm{\Omega }_c^2}},`$ (1a)
$`T_{fi}^{\left(3\right)}={\displaystyle \frac{\delta _{\sigma _L,X}\kappa _L\mathrm{\Omega }_c/\sqrt{2}}{\left(E_{1d}+\mathrm{}ck_LE_{3d}+i\mathrm{\Gamma }_3/2\right)\left(E_{1d}+\mathrm{}ck_LE_{2d}\right)\mathrm{\Omega }_c^2}}\times {\displaystyle \frac{\delta _{\sigma _R,X}\left|\kappa _R\right|^2/2}{E_{1d}+\mathrm{}ck_L^{}+\mathrm{}ck_R^{}E_{4d}+i\mathrm{\Gamma }_4/2}}`$
$`\times {\displaystyle \frac{\delta _{\sigma _L,X}\kappa _L^{}\mathrm{\Omega }_c^{}/\sqrt{2}}{\left(E_{1d}+\mathrm{}ck_L^{}E_{3d}+i\mathrm{\Gamma }_3/2\right)\left(E_{1d}+\mathrm{}ck_L^{}E_{2d}\right)\mathrm{\Omega }_c^2}},`$ (1b)
where $`\mathrm{\Gamma }_{3(4)}`$ is the decay rate of the polariton states $`|3`$ (or $`|4`$). In silicon microspheres, the whispering gallery modes can have Q-factor as high as $`10^8`$ Wang et al. (2000), so the intrinsic decay of the cavity modes can be neglected. The relaxation rate of trions is of the order of $`\mu `$eV, much less than the cavity-to-fiber loss. Thus, the decay of the trion polaritons is dominated by the leakage of the cavity modes into the fiber modes. The decay rates thus can be approximated as $`\mathrm{\Gamma }_{3(4)}\left|\kappa _{L(R)}\right|^2/c`$.

The $`Y`$-polarized photons are not scattered not being coupled to the cavities by design. The linear term in Eq. (1a) contributes to the reflection of the $`X`$-polarized photon from the left channel. It is significantly suppressed by the EIT effect Schmidt and Imamo$`\overline{g}`$lu (1996), which results from the destructive interference between the damped polariton state $`|3d`$ and the meta-stable state $`|2d`$ coupled by the classical pump field, as is evident from the vanishing of $`T_{fi}^{\left(1\right)}`$ when the incoming photon is in resonance with the $`|1d|3d`$ transition. The linear reflection of the $`X`$-polarized photons from the right channel is eliminated since the initial state of the system has been prepared in $`|1d`$ by the laser cooling cycle: The classical pulse pumps the ground state $`|x`$ to $`|3`$, and then the polariton state relax to $`|x+`$ through the cavity-to-fiber leakage (details to be published). The reduction of the linear reflection in both fibers brings to prominence the third order terms which are responsible for the gate operation. The nonlinear scattering term in Eq. (1b) is composed of three fractions corresponding to three processes: the excitation of the trion-polariton by the left-channel photon, the induced scattering of the right-channel photon, and the emission of the left-channel photon by the polariton recombination.

Due to the resonance features in the T-matrix, the transmission coefficient $`t_{XX}=fe^{i\varphi }`$, where $`f`$ and $`\varphi `$ are functions of the wavevectors, can cause the amplitude and phase modulation of the transmitted wave since the incoming photons are in wave packets. The amplitude modulation can be suppressed either by using longer time pulses or by working in far-resonance region. Although often overlooked in phase-shift estimation based on $`\chi ^{(3)}`$ susceptibility, the phase-variation effect results in distortion of the pulse shape and entanglement of the motion and polarization of the photons. The polarization states of the two photons are obtained by projection after the transmission. The effect of pulse deformation may be reduced by frequency filtering the transmitted pulse.

We show how shaping the input pulses leads to reduction of the output pulse deformation. (1) As a consequence of the optical coupling between the $`|3d`$ and $`|2d`$ states, the choice of the left input photon to be within $`\pm \mathrm{\Omega }_c`$ of being in resonance with the $`|1d|3d`$ transition, the linear reflection is reduced and the first factor on the right side of Eq. (1b) will yield a strong third-order transmission. (2) To diminish the pulse distortion due to the sharp resonant structure around the $`|2d|4d`$ transition, the right-channel pulse is detuned about $`\mathrm{\Gamma }_4/2`$ below the transition where the real part (corresponding to the phase-shift) of $`T^{(3)}`$ is large and flat while the imaginary part (corresponding to the reflection) has decreased to a small value. (3) To minimize the pulse broadening and distortion from the convolution of the input pulses with energy conservation, we choose the two input pulses to have sqaure-shaped spectra with much different widths. In our design, $`\mathrm{\Gamma }_4/2`$ is much larger than $`\mathrm{\Omega }_c`$, so the right-channel pulse is set the wider in frequency.

Figure 2 presents the transmitted wavefunction ($`k_L^{}`$, $`k_R^{}>0`$) for incoming photons in square pulses with the polarization state $`|XX`$. Though visible, the pulse distortion and broadening and the inhomogeneity in phase-shift is quite small. A conditional phase shift of $`\pi /29`$ is obtained with a transmission probability of $`0.72`$ and a 0.99 fidelity.

The loss in transmission and the pulse distortion in Fig. 2 result mainly from the imperfect EIT when the photon is off-resonant with the $`|1d|3d`$ transition. Improvement of both the pulse shape and transmission is effected by increasing the pump power $`\mathrm{\Omega }_c`$ (in order to open a larger EIT window) and by using narrower bandwidth pulses at the expense of a weak nonlinear phase shift. An example is shown in Fig. 3, in which a nonlinear phase shift $`\pi /330`$ is obtained almost without pulse-shape change or reflection loss, shown by the computed transmission probability of $`0.982`$ and almost perfect fidelity.

Although a small entanglement suffices for some quantum information purposes Cirac et al. (2001), the small phase-shift is not useful for most two-qubit operations. A large phase-shift can be achieved by using a series of many identical gates integrated into a single chip. With modern fabricating techniques, the integrated quantum gates can be constructed either with micro-disks and wave guides etched on semiconductor heterostructures Hagness et al. (1997) or with point- and line-defects engineered in photonic lattices Scherer et al. (2002).

To use the system to produce an entangled photon pair rather than to perform a controlled phase operation, we optimize the entanglement by a different procedure. The quantum operation is favored by maximizing the transmission but the entanglement is favored by symmetrizing the two photons for maximal projection of the polarization degrees of freedom. First the input state is prepared as the equal linear combination of the four polarization states of the two photons. Then the state is passed $`n`$ times through the coupled system as described above, undergoes single bit operation swapping the $`|X`$ and $`|Y`$ states in both photons, and is passed through the phase gate $`n`$ more times. This symmetrizes the transmitted density matrix after projecting out the motional degrees of freedom. Table 1 shows the calculated results for 2+2 and 4+4 gates. The transmission probabilities $`T`$ are much lower than for the phase operation. The quantitative measures of the operation including fidelity $`\mathrm{Tr}\left[\rho _t\rho _{ideal}\right]`$ towards the maximally entangled state $`\frac{1}{\sqrt{2}}(|XY+|YX)`$, the purity $`\mathrm{Tr}\left[\rho _t^2\right]`$, the concurrence $`C`$ and the entanglement of formation $`E(C)`$ Wootters (1998), all show excellent entanglement.

In summary, we have proposed a solid-state controlled phase gate for two photons. The flying qubits are conducted through fibers coupled to scattering centers composed of microcavities connected by a doped semiconductor nanodot. This allows a fiber implementation of quantum information processor. Calculated results show that the system is flexible as a phase gate as well as producing strong entanglement. The trions in doped nanodot used for nonlinear interaction here can be replaced by other electronic systems, such as biexcitons in an undoped nanodot (results will be published elsewhere), states in nanoclusters, or even some strong transitions in rare-earth impurities, e.g., the 4d-5f transition in $`Er^{2+}`$. The microcavity may be microspheres or defects in photonic lattices. The structure has unique features, such as small size, integrability, and stability, useful for quantum information and for scalable quantum computing.

This Work was supported by NSF DMR-0099572, ARDA/ARO DAAD19-02-1-0183, and QuIST/AFOSR F49620-01-1-0497. LJS thanks Y. Fainman for helpful discussions.


format: text

text:
# THE FOURTH MOMENT OF THE ZETA-FUNCTION

## 1. Introduction

A central place in analytic number theory is occupied by the Riemann zeta-function $`\zeta (s)`$, defined for $`\mathrm{}\mathrm{e}s>1`$ by

$$\zeta (s)=\underset{n=1}{\overset{\mathrm{}}{}}n^s=\underset{p\text{prime}}{}(1p^s)^1,$$

and otherwise by analytic continuation. It admits meromorphic continuation to the whole complex plane, its only singularity being the simple pole $`s=1`$ with residue 1. From the functional equation (see , )

$$\zeta (s)=\chi (s)\zeta (1s),\chi (s)=2^s\pi ^{s1}\mathrm{sin}\left(\frac{1}{2}\pi s\right)\mathrm{\Gamma }(1s),$$

which is valid for any complex $`s`$, it follows that $`\zeta (s)`$ has zeros at $`s=2,4,\mathrm{}`$ . These zeros are traditionally called the “trivial” zeros of $`\zeta (s)`$, to distinguish them from the complex zeros of $`\zeta (s)`$, of which the smallest ones (in absolute value) are $`\frac{1}{2}\pm 14.134725\mathrm{}i`$. It is well-known that all complex zeros of $`\zeta (s)`$ lie in the so-called “critical strip” $`0<\sigma =\mathrm{}\mathrm{e}s<1`$, and the Riemann Hypothesis (RH for short) is the conjecture, stated in 1859 by B. Riemann , that very likely all complex zeros of $`\zeta (s)`$ have real parts equal to $`\frac{1}{2}`$. For this reason the line $`\sigma =\frac{1}{2}`$ is called the “critical line” in the theory of $`\zeta (s)`$. The RH is extensively discussed by many authors, and recently by E. Bombieri .

The RH is now probably the most celebrated and difficult open problem in whole Mathematics. Its proof (or disproof) would have very important consequences in multiplicative number theory, especially in problems involving the distribution of primes.

## 2. Mean values of $`|\zeta (\frac{1}{2}+it)|`$

Mean values of $`|\zeta (\frac{1}{2}+it)|`$ are fundamental in the theory and applications of $`\zeta (s)`$ (see  for an extensive account). For $`k1`$ a fixed integer let

$$_0^T|\zeta (\frac{1}{2}+it)|^{2k}dt=TP_{k^2}(\mathrm{log}T)+E_k(T),$$
$`(2.1)`$
where for some suitable coefficients $`a_{j,k}`$ one has

$$P_{k^2}(y)=\underset{j=0}{\overset{k^2}{}}a_{j,k}y^j,$$
$`(2.2)`$
and in particular

$$P_1(y)=y+2\gamma 1\mathrm{log}(2\pi ),$$

where $`\gamma =\mathrm{\Gamma }^{}(1)=0.577\mathrm{}`$ is Euler’s constant. One hopes that

$$E_k(T)=o(T)(T\mathrm{})$$

will hold for every fixed integer $`k1`$, but so far this is known to be true only in the cases $`k=1`$ and $`k=2`$, when $`E_k(T)`$ is a true error term in (2.1). When $`k3`$ it is not even known what should be the values of the coeeficients $`a_{j,k}`$ in (2.2). The connection between $`E_k(T)`$ and the RH is indirect, namely there is a connection with the Lindelöf hypothesis (LH for short). The LH is also a famous unsettled problem, and it states that ($`fg`$ means $`|f(x)|<Cg(x)`$ for some $`C>0`$ and $`xx_0`$)

$$\zeta (\frac{1}{2}+it)_\epsilon t^\epsilon $$
$`(2.3)`$
for any given $`\epsilon >0`$ and $`tt_0>0`$ (since $`\overline{\zeta (\frac{1}{2}+it)}=\zeta (\frac{1}{2}it)`$, $`t`$ may be assumed to be positive). The RH implies LH, in fact it even gives (see ) an estimate stronger than (2.3), namely

$$\zeta (\frac{1}{2}+it)\mathrm{exp}\left(\frac{A\mathrm{log}t}{\mathrm{log}\mathrm{log}t}\right)(A>0,tt_0).$$

It is yet unknown whether the LH implies the RH. The best current unconditional bound for the order of $`\zeta (\frac{1}{2}+it)`$, which is far from the LH, is due to M.N. Huxley . This is

$$\zeta (\frac{1}{2}+it)_\epsilon t^{c+\epsilon },c=\frac{32}{205}=0.156098\mathrm{}.$$

Huxley  also proved

$$E_1(T)_\epsilon T^{\frac{137}{432}+\epsilon },\frac{137}{432}=0.31713\mathrm{}.$$

The LH is equivalent to the bound

$$_0^T|\zeta (\frac{1}{2}+it)|^{2k}dt_{k,\epsilon }T^{1+\epsilon }$$

for every $`k1`$ and any $`\epsilon >0`$, which in turn is the same as

$$E_k(T)_{k,\epsilon }T^{1+\epsilon }(kN).$$

We have $`\mathrm{\Omega }`$-results in the case $`k=1,2`$, which show that $`E_1(T)`$ and $`E_2(T)`$ cannot be always small. J.L. Hafner and the author  proved that

$$E_1(T)=\mathrm{\Omega }_+\left((T\mathrm{log}T)^{\frac{1}{4}}(\mathrm{log}\mathrm{log}T)^{\frac{3+\mathrm{log}4}{4}}\mathrm{e}^{C\sqrt{\mathrm{log}\mathrm{log}\mathrm{log}T}}\right)$$

and

$$E_1(T)=\mathrm{\Omega }_{}\left(T^{\frac{1}{4}}\mathrm{exp}\left(\frac{D(\mathrm{log}\mathrm{log}T)^{\frac{1}{4}}}{(\mathrm{log}\mathrm{log}\mathrm{log}T)^{\frac{3}{4}}}\right)\right)$$

for some absolute constants $`C,D>0`$. Moreover, the author  proved that there exist constants $`A,B>0`$ such that, for $`TT_0`$, every interval $`[T,T+B\sqrt{T}]`$ contains points $`t_1,t_2`$ for which

$$E_1(t_1)>At_1^{\frac{1}{4}},E_1(t_2)<At_2^{\frac{1}{4}}.$$

($`f=\mathrm{\Omega }_+(g)`$ means that $`lim\; supf/g>0`$, $`f=\mathrm{\Omega }_{}(g)`$ means that $`lim\; inff/g<0`$).

## 3. The fourth moment of $`|\zeta (\frac{1}{2}+it)|`$

The asymptotic formula for the fourth moment of $`\zeta (s)`$ on the critical line is customarily written as (this is (2.1) and (2.2) when $`k=2`$)

$$_0^T|\zeta (\frac{1}{2}+it)|^4dt=TP_4(\mathrm{log}T)+E_2(T),P_4(x)=\underset{j=0}{\overset{4}{}}a_jx^j.$$

A classical result of A.E. Ingham  from 1926 is that $`a_4=1/(2\pi ^2)`$ and that the error term $`E_2(T)`$ satisfies the bound $`E_2(T)T\mathrm{log}^3T`$. In 1979 D.R. Heath-Brown  made significant progress in this problem by proving that

$$E_2(T)_\epsilon T^{\frac{7}{8}+\epsilon }.$$

He also calculated

$$a_3=\mathrm{\hspace{0.33em}2}(4\gamma 1\mathrm{log}(2\pi )12\zeta ^{}(2)\pi ^2)\pi ^2$$

and produced more complicated expressions for $`a_0,a_1`$ and $`a_2`$. The author  made an explicit evaluation of $`a_0,a_1,a_2`$. For $`\mathrm{}\mathrm{e}s>1`$

$$\zeta ^2(s)=\underset{n=1}{\overset{\mathrm{}}{}}d(n)n^s,$$

where $`d(n)`$ is the sum of positive divisors of $`n`$. Note that $`\zeta ^4=\zeta ^2\zeta ^2`$ and that $`\zeta ^2`$ can be approximated by a finite sum $`_nd(n)n^s`$. Therefore by integrating

$$|\zeta (\frac{1}{2}+it)|^4=\zeta ^2(\frac{1}{2}+it)\zeta ^2(\frac{1}{2}it)$$

we are led to the asymptotic evaluation of the sum

$$\underset{nx}{}d(n)d(n+f),$$

where $`1fx`$ is not fixed. This is the so-called binary additive divisor problem (see ,  and ). Modern approaches both to the study of $`E_2(T)`$ and the binary additive divisor problem involve the use of spectral theory.

## 4. Spectral theory

For a competent account on spectral theory and its applications to $`\zeta (s)`$ the reader is referred to Y. Motohashi’s monograph . Here we shall only briefly present some basic facts. On the upper complex half-plane $`H`$ the modular group

$$\mathrm{\Gamma }=\mathrm{S}L(2,Z)=\{\left(\genfrac{}{}{0pt}{}{ab}{cd}\right):(a,b,c,dZ)(adbc=1)\}$$

acts by $`\gamma z=(az+b)/(cz+d)`$ if $`\gamma \mathrm{\Gamma }`$. The non-Euclidean Laplace operator

$$L=y^2\left(\frac{^2}{x^2}+\frac{^2}{y^2}\right)$$

is invariant under $`\mathrm{\Gamma }`$, as is the measure $`\mathrm{d}\mu (z)=y^2\mathrm{d}x\mathrm{d}y(z=x+iy).`$ Non-holomorphic cusp forms (so-called Maass wave forms) are eigenfunctions $`\mathrm{\Psi }(z)`$ of the discrete spectrum of $`L`$, which has the form

$$\left\{\lambda _j\right\}_{j=1}^{\mathrm{}},\lambda _j=\kappa _j^2+\frac{1}{4}(\kappa _j>0).$$

They satisfy $`L\mathrm{\Psi }(z)=\lambda \mathrm{\Psi }(z)`$, $`\mathrm{\Psi }(\gamma z)=\mathrm{\Psi }(z)`$ for $`\gamma \mathrm{\Gamma }`$, and the finiteness condition

$$_𝒟|\mathrm{\Psi }(z)|^2d\mu (z)=_𝒟|\mathrm{\Psi }(x+iy)|^2y^2dxdy<\mathrm{}.$$

Here $`𝒟`$ is the fundamental domain of $`\mathrm{\Gamma }`$, namely

$$𝒟=\{z:y>0,|z|>1,\frac{1}{2}x\frac{1}{2}\}\{z:|z|=1,\frac{1}{2}x0\}.$$

Let henceforth $`\phi _j`$ be the Maass wave form attached to $`\kappa _j`$ so that $`\{\phi _j\}_{j=1}^{\mathrm{}}`$ forms an orthonormal basis with respect to the Petersson inner product

$$(f_1,f_2):=_𝒟f_1\overline{f}_2d\mu (z),$$

and $`\phi _j`$ is an eigenfunction of every Hecke operator. The Hecke operator $`T_n`$ acts on $`H`$, for given $`nN`$, by the relation

$$(T_nf)(z)=n^{1/2}\underset{ad=n,d>0}{}\underset{b(\text{mod}d)}{}f\left(\frac{az+b}{d}\right).$$

We have the Fourier expansion

$$\phi _j(z)=\underset{n0}{}\rho _j(n)e^{2\pi ix}\sqrt{y}K_{i\kappa _j}(2\pi |n|y)$$

for $`z=x+iy,\rho _j(n)=\overline{\rho _j(n)},`$ where $`K`$ is the Bessel function (the Mcdonald function)

$$K_s(z)=\frac{1}{2}_0^{\mathrm{}}t^{s1}\mathrm{exp}\left(\frac{1}{2}z(t+\frac{1}{t})\right)dt(\mathrm{}\mathrm{e}z>0).$$

For $`nN`$, $`t_j(n)`$ is the eigenvalue corresponding to $`\phi _j`$ with respect to $`T_n`$,

$$T_n\phi _j(z)=t_j(n)\phi _j(z).$$

Then

$$\rho _j(1)t_j(n)=\rho _j(n).$$

The Hecke series attached to $`\phi _j(z)`$ is

$$H_j(s):=\underset{n=1}{\overset{\mathrm{}}{}}t_j(n)n^s=\underset{p}{}(1t_j(p)p^s+p^{2s})^1(\mathrm{}\mathrm{e}s>1).$$

It is known that $`H_j(s)`$ continues analytically to an entire function over $`C`$, and that for any $`s`$ satisfies the functional equation

$$H_j(s)=\pi ^1(2\pi )^{2s1}\mathrm{\Gamma }(1s+i\kappa _j)\mathrm{\Gamma }(1si\kappa _j)\left\{\mathrm{cos}(\pi s)+\epsilon _j\mathrm{cosh}(\pi \kappa _j)\right\}H_j(1s)$$

with $`\epsilon _j=\pm 1`$.

## 5. Motohashi’s explicit formula

Recent progress on the fourth moment of $`|\zeta (\frac{1}{2}+it)|`$ is primarily due to the fundamental explicit formula of Y. Motohashi (see , , , ) of 1989 for

$$I(T,\mathrm{\Delta }):=\frac{1}{\mathrm{\Delta }\sqrt{\pi }}_{\mathrm{}}^{\mathrm{}}|\zeta (\frac{1}{2}+iT+it)|^4\mathrm{e}^{t^2\mathrm{\Delta }^2}dt,$$

under the condition

$$0<\mathrm{\Delta }\frac{T}{\mathrm{log}T}.$$

The presence of the Gaussian smoothing factor $`\mathrm{e}^{x^2}`$ enabled Motohashi to deal with various convergence problems occurring in the proof. In simplified form, the formula is

$$\begin{array}{cc}\hfill I(T,\mathrm{\Delta })& =\pi 2^{1/2}T^{1/2}\underset{j=1}{\overset{\mathrm{}}{}}\alpha _j\kappa _j^{1/2}H_j^3(\frac{1}{2})\mathrm{sin}\left(\kappa _j\mathrm{log}\left(\frac{\kappa _j}{4\mathrm{e}T}\right)\right)\mathrm{e}^{(\mathrm{\Delta }\kappa _j/2T)^2}\hfill \\ & +O(\mathrm{log}^CT).\hfill \end{array}$$
$`(5.1)`$
Here

$$\alpha _j=\frac{|\rho _j(1)|^2}{\mathrm{cosh}(\pi \kappa _j)},$$
$`(5.2)`$
and, for arbitrary fixed $`A>0`$, $`C=C(A)>0`$ with $`\mathrm{\Delta }`$ satisfying

$$\frac{\sqrt{T}}{\mathrm{log}^AT}\mathrm{\Delta }T\mathrm{exp}(\sqrt{\mathrm{log}T}).$$

The proof of the formula for $`I(T,\mathrm{\Delta })`$ depends heavily on the use of spectral theory, especially the so-called Kuznetsov trace formulas, which relate sums of Kloosterman sums (see (6.4)) to certain sums involving quantities from spectral theory.

## 6. New results on the fourth moment

At present the bound

$$_0^T|\zeta (\frac{1}{2}+it)|^kdt_\epsilon T^{1+\epsilon }$$

is not known to hold for any constant $`k>4`$. Therefore the function $`E_2(T)`$ is particularly important in the theory of mean values of $`\zeta (s)`$. In recent years much advance has been made in connection with $`E_2(T)`$ and related problems. One of the main problems is to get rid of the Gaussian smoothing factor $`\mathrm{e}^{x^2}`$ in Motohashi’s formula (5.1) and apply it to obtain results on $`E_2(T)`$ itself.

Motohashi and the author in four papers – obtained several results on $`E_2(T)`$ and some related problems. It is known now that ($`f=\mathrm{\Omega }(g)`$ means that $`lim\; sup|f|/g>0`$)

$$E_2(T)=O(T^{2/3}\mathrm{log}^{C_1}T),E_2(T)=\mathrm{\Omega }(T^{1/2}),$$

$$\underset{0}{\overset{T}{}}E_2(t)dt=O(T^{3/2}),\underset{0}{\overset{T}{}}E_2^2(t)dt=O(T^2\mathrm{log}^{C_2}T),$$

with effective constants $`C_1,C_2>0`$ (the values $`C_1=8,C_2=22`$ are admissible). Y. Motohashi  improved the omega-result to

$$E_2(T)=\mathrm{\Omega }_\pm (T^{1/2}).$$

Finally the author  made further progress in this problem by proving the following quantitative omega-result: there exist two constants $`A>0,B>1`$ such that for $`TT_0>0`$ every interval $`[T,BT]`$ contains points $`T_1,T_2`$ for which

$$E_2(T_1)>AT_1^{1/2},E_2(T_2)<AT_2^{1/2}.$$

This follows from the asymptotic formula

$$\begin{array}{cc}\hfill \underset{0}{\overset{\mathrm{}}{}}E_2(t)\mathrm{e}^{t/T}dt& =2T^{\frac{3}{2}}\mathrm{}\mathrm{e}\left\{\underset{j=1}{\overset{\mathrm{}}{}}\alpha _jH_j^3(\frac{1}{2})R(\kappa _j)\mathrm{\Gamma }(\frac{1}{2}i\kappa _j)T^{\kappa _j}\right\}\hfill \\ & +O\left(T^{\frac{3}{2}}\mathrm{exp}\left(\frac{\mathrm{log}T}{\mathrm{log}\mathrm{log}T}\right)\right),\hfill \end{array}$$
$`(6.1)`$
with $`\alpha _j`$ given by (5.2) and

$$R(y):=\sqrt{\frac{\pi }{2}}\left(2^{iy}\frac{\mathrm{\Gamma }(\frac{1}{4}+\frac{i}{2}y)}{\mathrm{\Gamma }(\frac{1}{4}\frac{i}{2}y)}\right)^3\mathrm{\Gamma }(2iy)\mathrm{cosh}(\pi y).$$
$`(6.2)`$
Y. Motohashi and the author  proved that

$$\underset{r=1}{\overset{R}{}}_{t_r}^{t_r+\mathrm{\Delta }}|\zeta (\frac{1}{2}+it)|^4dtR\mathrm{\Delta }\mathrm{log}^4T+TR^{1/2}\mathrm{\Delta }^{1/2}\mathrm{log}^CT$$

for some $`C>0`$, where $`\mathrm{log}T\mathrm{\Delta }\frac{T}{\mathrm{log}T}`$ and

$$Tt_1<t_2<\mathrm{}<t_R2T,t_{r+1}t_r\mathrm{\Delta }(r=1,\mathrm{},R1).$$

From this result it is not difficult to obtain as a corollary the bound

$$_0^T|\zeta (\frac{1}{2}+it)|^{12}dt=O(T^2\mathrm{log}^BT),$$

proved first (with $`B=17`$) by Heath-Brown  in 1978.

The asymptotic formula (6.1) for $`_0^{\mathrm{}}E_2(t)\mathrm{e}^{t/T}dt`$ is a Laplace transform formula. Integral transforms such as the Mellin transform (see , ) and the Laplace transform (see , ) play an important rôle in analytic number theory, in particular in the theory of $`\zeta (s)`$. One can consider the general function

$$L_k(s):=\underset{0}{\overset{\mathrm{}}{}}|\zeta (\frac{1}{2}+ix)|^{2k}e^{sx}dx(kN,\mathrm{}\mathrm{e}s>0).$$

A classical result of H. Kober  from 1936 says that, as $`\sigma 0+`$,

$$L_1(2\sigma )=\frac{\gamma \mathrm{log}(4\pi \sigma )}{2\mathrm{sin}\sigma }+\underset{n=0}{\overset{N}{}}c_n\sigma ^n+O(\sigma ^{N+1})$$

for any given integer $`N1`$, where the $`c_n`$’s are effectively computable constants.

F.V. Atkinson  obtained in 1941 the asymptotic formula

$$L_2(\sigma )=\frac{1}{\sigma }\left(A\mathrm{log}^4\frac{1}{\sigma }+B\mathrm{log}^3\frac{1}{\sigma }+C\mathrm{log}^2\frac{1}{\sigma }+D\mathrm{log}\frac{1}{\sigma }+E\right)+\lambda _2(\sigma ),$$
$`(6.3)`$
where $`\sigma 0+`$,

$$A=\frac{1}{2\pi ^2},B=\pi ^2(2\mathrm{log}(2\pi )6\gamma +24\zeta ^{}(2)\pi ^2),\lambda _2(\sigma )_\epsilon \left(\frac{1}{\sigma }\right)^{\frac{13}{14}+\epsilon }.$$

He also indicated how, by the use of estimates for Kloosterman sums

$$S(m,n;c):=\underset{1d<c,(d,c)=1,dd^{}1(\text{mod}c)}{}\mathrm{e}\left(\frac{md+nd^{}}{c}\right)(\mathrm{e}(z)=\mathrm{e}^{2\pi iz}),$$
$`(6.4)`$
one can improve the exponent $`\frac{13}{14}`$ to $`\frac{8}{9}`$. This is important historically, in view of contemporary importance of Kloosterman sums which stems from the Kuznetsov trace formulas and other important applications.

The author  gave explicit, albeit complicated expressions for the remaining coefficients $`C,D`$ and $`E`$ in (6.3) and improved Atkinson’s bound for $`\lambda _2(\sigma )`$ to

$$\lambda _2(\sigma )\sigma ^{1/2}(\sigma 0+).$$

Recently in  he proved a result which generalizes and sharpens (6.3): Let $`0\varphi <\frac{\pi }{2}`$ be given. Then for $`0<|s|1`$ and $`|\mathrm{arg}s|\varphi `$ we have

$$\begin{array}{cc}\hfill L_2(s)& =\frac{1}{s}(A\mathrm{log}^4\frac{1}{s}+B\mathrm{log}^3\frac{1}{s}+C\mathrm{log}^2\frac{1}{s}+D\mathrm{log}\frac{1}{s}+E)\hfill \\ & +s^{\frac{1}{2}}\left\{\underset{j=1}{\overset{\mathrm{}}{}}\alpha _jH_j^3(\frac{1}{2})\left(s^{i\kappa _j}R(\kappa _j)\mathrm{\Gamma }(\kappa _j)+s^{i\kappa _j}R(\kappa _j)\mathrm{\Gamma }(\kappa _j)\right)\right\}+G_2(s),\hfill \end{array}$$

where the constants $`A,\mathrm{},E`$ are as in (6.3), $`R(y)`$ is given by (6.2) and in the above region $`G_2(s)`$ is a regular function satisfying ($`C>0`$ is a suitable constant)

$$G_2(s)|s|^{\frac{1}{2}}\mathrm{exp}\left\{\frac{C\mathrm{log}(|s|^1)}{(\mathrm{log}\mathrm{log}(|s|^1))^{\frac{2}{3}}(\mathrm{log}\mathrm{log}\mathrm{log}(|s|^1))^{\frac{1}{3}}}\right\}.$$

In  the author proved: There exist constants $`A>0`$ and $`B>1`$ such that, for $`TT_0>0`$, every interval $`[T,BT]`$ contains points $`t_1,t_2`$ for which

$$_0^{t_1}E_2(t)dt>At_1^{3/2},_0^{t_2}E_2(t)dt<At_2^{3/2}.$$

This result implies that

$$_0^TE_2(t)dt=\mathrm{\Omega }_\pm (T^{3/2}),$$

and that

$$\underset{n\mathrm{}}{lim\; sup}\frac{\mathrm{log}(u_{n+1}u_n)}{\mathrm{log}u_n}\mathrm{\hspace{0.33em}1},$$

where $`u_n`$ is the $`n`$–th zero of $`E_2(T)`$. It can also be used to prove the lower bound result

$$_0^TE_2^2(t)dtT^2,$$

which complements the earlier result with Motohashi (cf. ) that

$$_0^TE_2^2(t)dtT^2\mathrm{log}^CT(C=22).$$

Another recent result of the author , which provides an asymptotic formula for the integral of $`E_2(t)`$, is

$$\begin{array}{cc}& _0^TE_2(t)dt=O\{T^{\frac{3}{2}}\mathrm{exp}(C(\mathrm{log}T)^{\frac{3}{5}}(\mathrm{log}\mathrm{log}T)^{\frac{1}{5}})\}+\hfill \\ & +\mathrm{\hspace{0.17em}2}T^{\frac{3}{2}}\mathrm{}\mathrm{e}\left\{\underset{j=1}{\overset{\mathrm{}}{}}\alpha _jH_j^3(\frac{1}{2})\frac{T^{i\kappa _j}}{(\frac{1}{2}+i\kappa _j)(\frac{3}{2}+i\kappa _j)}R(\kappa _j)\right\}.\hfill \end{array}$$

Here $`C>0`$ is a suitable constant, and the error term depends on the best known zero-free region for $`\zeta (s)`$ (see \[9, Chapter 6\]).

In spite of significant recent results on $`E_2(T)`$, many open problems remain. Here are four of them.

Problem 1. Does there exist $`A>0`$ such that, as $`T\mathrm{}`$,

$$_0^TE_2^2(t)dtAT^2\mathrm{?}$$

Problem 2. Is it true that for every $`\epsilon >0`$

$$E_2(T)=O_\epsilon \left(T^{\frac{1}{2}+\epsilon }\right)\mathrm{?}$$

Problem 3. Does one have

$$\underset{T\mathrm{}}{lim\; sup}|E_2(T)|T^{\frac{1}{2}}=\mathrm{}\mathrm{?}$$

Problem 4. Does one have

$$\underset{n\mathrm{}}{lim\; sup}\frac{\mathrm{log}(u_{n+1}u_n)}{\mathrm{log}u_n}<\mathrm{\hspace{0.33em}1}\mathrm{?}$$

2cm

 F.V. Atkinson, ‘The mean value of the zeta-function on the critical line’, Quart. J. Math. Oxford 10(1939), 122-128.

  F.V. Atkinson, ‘The mean value of the zeta-function on the critical line’, Proc. London Math. Soc. (2)47(1941), 174-200.

 E. Bombieri, ‘Problems of the Millenium: the Riemann Hypothesis’, 2000, http//www.ams.org /claymath/prize\_problems/riemann.pdf., 11pp.

 J.L. Hafner and A. Ivić, ‘On the mean square of the Riemann zeta-function on the critical line’, J. Number Theory 32(1989), 151-191.

 D.R. Heath-Brown, ‘The twelfth power moment of the Riemann zeta-function’, Quart. J. Math. Oxford 29(1978), 443-462.

  D.R. Heath-Brown, ‘The fourth moment of the Riemann zeta-function’, Proc. London Math. Soc. (3)38(1979), 385-422.

 M.N. Huxley, ‘Integer points, exponential sums and the Riemann zeta-function’, Proc. of “Millenial Conference on Number Theory”, Urbana, May 2000, in print.

  A.E. Ingham, ‘Mean-value theorems in the theory of the Riemann zeta-function’, Proc. London Math. Soc. (2)27(1926), 273-300.

  A. Ivić, ‘The Riemann zeta-function’, John Wiley and Sons, New York, 1985.

 A. Ivić, ‘Large values of certain number-theoretic error terms’, Acta Arith. 56 (1990), 135-159.

  A. Ivić, ‘Mean values of the Riemann zeta-function’, LN’s 82, Tata Institute of Fundamental Research, Bombay, 1991 (distr. by Springer Verlag, Berlin etc.).

  A. Ivić, ‘On the fourth moment of the Riemann zeta-function’, Publs. Inst. Math. (Belgrade) 57(71)(1995), 101-110.

  A. Ivić, ‘The Mellin transform and the Riemann zeta-function’, in “Proceedings of the Conference on Elementary and Analytic Number Theory (Vienna, July 18-20, 1996)”, Universität Wien & Universität für Bodenkultur, Eds. W.G. Nowak and J. Schoißengeier, Vienna 1996, 112-127.

 A. Ivić, ‘On the error term for the fourth moment of the Riemann zeta-function’, J. London Math. Soc., 60(2)(1999), 21-32.

 A. Ivić, ‘On the integral of the error term in the fourth moment of the Riemann zeta-function’, Functiones et Approximatio, 28(2000), 37-48.

 A. Ivić, ‘The Laplace transform of the fourth moment of the zeta-function’, Univ. Beog. Publik. Elektroteh. Fak. Ser. Matematika 11(2000), 41-48.

 A. Ivić and Y. Motohashi, ‘A note on the mean value of the zeta and L-functions VII’, Proc. Japan Acad. Ser. A 66(1990), 150-152.

  A. Ivić and Y. Motohashi, ‘The mean square of the error term for the fourth moment of the zeta-function’, Proc. London Math. Soc. (3)66(1994), 309-329.

  A. Ivić and Y. Motohashi, ‘The fourth moment of the Riemann zeta-function’, J. Number Theory 51(1995), 16-45.

  A. Ivić and Y. Motohashi, ‘On some estimates involving the binary additive problem’, Quart. J. Math. Oxford (2)46(1995), 471-483.

 A. Ivić, M. Jutila and Y. Motohashi, The Mellin transform of powers of the zeta-function, Acta Arith. 95(2000), 305-342.

  M. Jutila, ‘The fourth moment of Riemann’s zeta-function and the additive divisor problem’, in “Analytic Number Theory; Proc. Conf. in Honor of Heini Halberstam, Vol. 2” (eds. B.C. Berndt et al.), Birkhäuser, Boston etc., 1996, 517-536.

  M. Jutila, ‘Mean values of Dirichlet series via Laplace transforms’, in “Analytic Number Theory” (ed. Y. Motohashi), London Math. Soc. LNS 247, Cambridge University Press, Cambridge, 1997, 169-207.

 H. Kober, Eine Mittelwertformel der Riemannschen Zetafunktion, Compositio Math. 3(1936), 174-189.

 Y. Motohashi, The fourth power mean of the Riemann zeta-function, in Proceedings of the Amalfi Conference on Analytic Number Theory 1989, eds. E. Bombieri et al., Università di Salerno, Salerno, 1992, 325-344.

 Y. Motohashi, ‘An explicit formula for the fourth power mean of the Riemann zeta-function’, Acta Math. 170(1993), 181-220.

 Y. Motohashi,‘The binary additive divisor problem’, Ann. Sci. École. Norm. Sup. $`4^e`$ série 27(1994), 529-572.

 Y. Motohashi, ‘A relation between the Riemann zeta-function and the hyperbolic Laplacian’, Annali Scuola Norm. Sup. Pisa, Cl. Sci. IV ser. 22(1995), 299-313.

 Y. Motohashi, ‘The Riemann zeta-function and the non-Euclidean Laplacian’, Sugaku Expositions, AMS 8(1995), 59-87.

 Y. Motohashi, ‘Spectral theory of the Riemann zeta-function’, Cambridge University Press, Cambridge, 1997.

 B. Riemann, ‘Über die Anzahl der Primzahlen unter einer gegebener Grösse’, Monatsber. Akad. Berlin (1859), 671-680.

 E.C. Titchmarsh, The Theory of the Riemann Zeta-Function, Clarendon Press, Oxford, 1951.

Aleksandar Ivić

Katedra Matematike RGF-a

Universitet u Beogradu, -Dušina 7

11000 Beograd, Serbia (Yugoslavia)

aivicrgf.bg.ac.yu, aivicmatf.bg.ac.yu


text:
# Rapid LISA Astronomy

## Abstract

A simple method is presented for removing the amplitude, frequency and phase modulations from the Laser Interferometer Space Antenna (LISA) data stream for sources at any sky location. When combined with an excess power trigger or the fast chirp transform, the total demodulation procedure allows the majority of LISA sources to be identified without recourse to matched filtering.

A revolution is underway that will transform the field of astronomy. An array of ground and space based gravitational wave detectors are in various stages of design, construction and operation. The first space based detector, the Laser Interferometer Space Antenna lppa , is being developed for a launch in 2011. These gravitational wave detectors will open a new window on the Universe that complements traditional forms of astronomy hughes .

However, before we can talk about gravitational wave astronomy, a method must be found to turn a gravitational wave detector into a gravitational wave telescope - that is, a device that is capable of locating individual sources on the sky and determining their physical characteristics. In traditional electromagnetic astronomy the conversion of a detector - such as the human eye or a photographic plate - into an observational instrument is achieved by building a device that collects and focuses the incident radiation onto the detector. Such an approach is impractical for gravitational radiation as the weak coupling of the radiation to matter prevents the construction of a physical telescope. Instead, a gravitational wave detector is transformed into a gravitational wave telescope by data analysis algorithms.

In principle we would like to use matched filtering of the data against suitable waveform templates as this delivers the highest fidelity results, but in practice the computational cost of matched filtering can be prohibitive. The solution is to adopt a hierarchical approach where a fast, but sub-optimal, method is used to find a good initial solution, which can then be refined using matched filtering. The excess power statistic lads  used in ground based gravitational wave astronomy is a good example of this approach. However, we cannot simply take the ground based data analysis algorithms and apply them to the space based systems as the high frequency, ground based detectors primarily look for short lived, burst like signals, whereas the low frequency space based detectors primarily look for long lived, continuous sources. One of the key differences between ground and space based gravitational wave detection is ratio of the observation time to the orbital period of the detector. The ratio is small for most high frequency sources, and of order one for most low frequency sources. This means that LISA data analysis algorithms have to deal with orbital signal modulation and multiple overlapping sources.

Here we describe a simple method for converting the output of a space based gravitational wave detector, such as LISA, into useful astronomical information about the location and physical characteristics a gravitational wave source. The method works by removing the amplitude, phase and frequency modulation imparted by the orbital motion of the detector. For monochromatic sources, the demodulation procedure focuses the signal into a narrow frequency band, which creates a spike in the demodulated power spectrum. By searching for sky locations that yield the largest spikes it is possible to locate the brightest sources on the sky. These can then be subtracted and the procedure iterated until all resolvable sources have been found. For chirping sources the demodulation can be followed by a fast chirp transform tom , which effectively re-concentrates the power that is spread by the frequency evolution of the source. Once again the brightest sources on the sky can be identified and iteratively subtracted to expose the weaker sources.

The demodulation procedure is remarkably simple. It utilizes the fact that space based detectors such as LISA are able to return several interferometer outputs aet . A general gravitational wave coming from a source in the $`\widehat{n}`$ direction can be expressed in barycentric coordinates as

$$h_{\mu \nu }(t,\widehat{n})=h_+(t,\widehat{n})e_{\mu \nu }^+(\widehat{n})+h_\times (t,\widehat{n})e_{\mu \nu }^\times (\widehat{n})$$
(1)
where $`h_+(t,\widehat{n})`$ and $`h_\times (t,\widehat{n})`$ describe the two polarization states with basis tensors $`e_{\mu \nu }^+(\widehat{n})`$ and $`e_{\mu \nu }^\times (\widehat{n})`$. The polarization states $`h_+`$, $`h_\times `$ are related to the principal polarization states $`h_+^P`$, $`h_\times ^P`$ by the polarization angle $`\psi `$:

$`h_+=h_+^P\mathrm{cos}2\psi +h_\times ^P\mathrm{sin}2\psi `$
$`h_\times =h_\times ^P\mathrm{cos}2\psi h_+^P\mathrm{sin}2\psi .`$ (2)
In the idealized limit of a single low frequency source and a noise-free detector, the interferometer outputs $`s(t)`$ can be expressed as

$$s(t)=\overline{h}_+(t,\widehat{n})D^+(t,\widehat{n})+\overline{h}_\times (t,\widehat{n})D^\times (t,\widehat{n})$$
(3)
where $`D^+(t,\widehat{n})`$ and $`D^\times (t,\widehat{n})`$ are the sky-location dependent interferometer response functions cr1  and $`\overline{h}_+(t,\widehat{n})`$ and $`\overline{h}_\times (t,\widehat{n})`$ are the Doppler shifted gravitational wave strains

$$\overline{h}(t,\widehat{n})=h(t+\frac{R}{c}\mathrm{sin}\theta \mathrm{cos}(2\pi f_mt\varphi ),\widehat{n}).$$
(4)
Here $`\widehat{n}(\theta ,\varphi )`$ denotes the location of the source in ecliptic coordinates, $`R`$ is the Earth-Sun distance and $`f_m=1/\mathrm{year}`$ is the modulation frequency. A simple, but impractical, approach for demodulating the signal is to form the combinations

$$\overline{h}_\times (t,\widehat{n})=\frac{s_\mathrm{I}(t)D_{\mathrm{II}}^+(t,\widehat{n})s_{\mathrm{II}}(t)D_\mathrm{I}^+(t,\widehat{n})}{D_\mathrm{I}^\times (t,\widehat{n})D_{\mathrm{II}}^+(t,\widehat{n})D_{\mathrm{II}}^\times (t,\widehat{n})D_\mathrm{I}^+(t,\widehat{n})}$$
(5)
and

$$\overline{h}_+(t,\widehat{n})=\frac{s_\mathrm{I}(t)D_{\mathrm{II}}^\times (t,\widehat{n})s_{\mathrm{II}}(t)D_\mathrm{I}^\times (t,\widehat{n})}{D_\mathrm{I}^\times (t,\widehat{n})D_{\mathrm{II}}^+(t,\widehat{n})D_{\mathrm{II}}^\times (t,\widehat{n})D_\mathrm{I}^+(t,\widehat{n})}$$
(6)
from interferometer outputs $`\mathrm{I}`$ and $`\mathrm{II}`$, then perform the coordinate transformation

$$t=t^{}\frac{R}{c}\mathrm{sin}\theta \mathrm{cos}(2\pi f_mt^{}\varphi )$$
(7)
to remove the leading order Doppler modulation cl1 , thereby recovering the demodulated strains $`\overline{h}_+(t^{})`$ and $`\overline{h}_\times (t^{})`$. The main problem with this approach is that the denominators in Eqns. (5) and (6) vanish at certain times of the year for sources within $`60`$ degrees of the ecliptic plane. When detector noise or other sources are present the inversion becomes singular. Although this simple demodulation procedure does not work in the time domain, it can be made to work in the frequency domain.

The key advantage of moving to the frequency domain is that the singularities in the temporal inversion do not correspond to singularities in the frequency space inversion. By Fourier transforming the data over a finite time interval, the demodulation procedure can be re-cast as a linear algebra problem that can be solved using the robust technique of singular value decomposition.

To simply our expressions we assume that the observation time $`T`$ is equal to one year so that the frequency resolution $`\mathrm{\Delta }f=1/T`$ is equal to the modulation frequency $`f_m`$ (the general expressions are not hard to derive, but they are more complicated). The Fourier space equivalents of (3) and (4) become

$$s_k=A_{kl}^+(\widehat{n})\overline{h}_l^+(\widehat{n})+A_{kl}^\times (\widehat{n})\overline{h}_l^\times (\widehat{n}),$$
(8)
and

$$\overline{h}_l^+=B_{ln}(\widehat{n})h_n^+(\widehat{n}),\overline{h}_l^\times =B_{ln}(\widehat{n})h_n^\times (\widehat{n}),$$
(9)
respectively. The matrices $`𝐀`$ and $`𝐁`$ are given by

$$A_{kl}^+(\widehat{n})=D_{kl}^+(\widehat{n}),A_{kl}^\times (\widehat{n})=D_{kl}^\times (\widehat{n})$$
(10)
and

$$B_{ln}(\widehat{n})=J_{(ln)}\left(2\pi nf_m\frac{R}{c}\mathrm{sin}\theta \right)e^{i(ln)(\pi /2\varphi )}.$$
(11)
Here $`J_q(x)`$ is a Bessel function of the first kind of order $`q`$ and we have used the Einstein summation convention for repeated indices. It is important to note that the matrices $`𝐀`$ and $`𝐁`$ are both band diagonal as $`D_j^+(\widehat{n})0D_j^\times (\widehat{n})`$ for $`|j|>4`$ and $`J_q(x)0`$ for $`|q|>x`$. The next step is to interleave the Fourier coefficients of the two polarizations $`h^+(t)`$ and $`h^\times (t)`$ into a single column vector $`𝐡`$ such that $`h_{2n}=h_n^+`$ and $`h_{2n+1}=h_n^\times `$. Similarly, the Fourier coefficients of the two interferometer outputs $`s_\mathrm{I}(t)`$ and $`s_{\mathrm{II}}(t)`$ are interleaved into a single column vector $`𝐬`$. The detector response can then be written as

$$𝐬=𝐌(\widehat{n})𝐡(\widehat{n}),$$
(12)
where the modulation matrix $`𝐌(\widehat{n})`$ is the interleaved product of $`𝐀^+(\widehat{n})`$, $`𝐀^\times (\widehat{n})`$ and $`𝐁(\widehat{n})`$:

$`M_{2k\mathrm{\hspace{0.17em}2}l}=A_{\mathrm{I}}^{+}{}_{kj}{}^{}B_{jl},M_{2k\mathrm{\hspace{0.17em}2}l+1}=A_{\mathrm{I}}^{\times }{}_{kj}{}^{}B_{jl},`$
$`M_{2k+\mathrm{1\hspace{0.17em}2}l}=A_{\mathrm{II}}^{+}{}_{kj}{}^{}B_{jl},M_{2k+\mathrm{1\hspace{0.17em}2}l+1}=A_{\mathrm{II}}^{\times }{}_{kj}{}^{}B_{jl}.`$ (13)
If the interferometer outputs are sampled $`N`$ times in a year, $`𝐌(\widehat{n})`$ will be a $`2N\times 2N`$ dimensional, complex, band diagonal matrix. Multiplying $`𝐬`$ by the inverse of the modulation matrix demodulates any sources located in the $`\widehat{n}`$ direction. In practice, $`𝐬`$ will contain contributions from instrument noise $`𝐧`$, and multiple sources $`𝐡^i(\widehat{n}^i)`$:

$$𝐬=\underset{i}{}𝐌^i(\widehat{n}^i)𝐡^i(\widehat{n}^i)+𝐧.$$
(14)
We can estimate the contribution to the signal $`𝐬`$ from a source located in the $`\widehat{n}`$ direction by solving for the vector $`𝐡^{\mathrm{eff}}(\widehat{n})`$ that satisfies the equation

$$𝐬=𝐌(\widehat{n})𝐡^{\mathrm{eff}}(\widehat{n}).$$
(15)
The best fit solution is found using a singular value decomposition. Suppose that source $`j`$ happens to lie in the $`\widehat{n}`$ direction. Then

$$𝐡^{\mathrm{eff}}(\widehat{n})𝐡^j(\widehat{n})+\underset{ij}{}𝐌^1(\widehat{n})𝐌^i(\widehat{n}^i)𝐡^i(\widehat{n}^i)+𝐌^1(\widehat{n})𝐧.$$
(16)
The key to the method is that sources at other sky locations are not demodulated, and the noise is merely reshuffled, not amplified. Thus, we can use the solution vector $`𝐡^{\mathrm{eff}}(\widehat{n})`$ as an estimate for the source vector $`𝐡^j(\widehat{n})`$. Another key point is that the modulation matrix is band diagonal, which allows us to perform the demodulation across limited bandwidths using small modulation sub-matrices. If this were not the case we would be confronted with the problem of inverting a $`10^8\times 10^8`$ matrix. Moreover, a source only has to compete with other sources that overlap with it in frequency space.

The demodulation procedure works for any type of gravitational wave source, but it is especially useful when the source is a mildly eccentric, nearly Newtonian binary as the demodulated signal is then approximately monochromatic. We can search for these sources by looking for spikes in the the power spectrum of $`𝐡^{\mathrm{eff}}(\widehat{n})`$ for different sky locations $`\widehat{n}`$. A similar approach was employed in Ref. cl1  using Doppler demodulation of a single interferometer output. Since the vast majority of LISA sources are non-chirping, circular Newtonian binaries, the combined demodulation, power spike search goes a long way toward solving the LISA data analysis problem. After establishing the sky location and frequency of a source, the next step is to refine the frequency measurement and extract the Fourier amplitudes of the two polarizations $`\stackrel{~}{h}^+`$ and $`\stackrel{~}{h}^\times `$. Unless the source frequency $`f`$ is an integer multiple of the sample frequency $`\mathrm{\Delta }f=f_m`$, the power will be shared by several discrete Fourier coefficients:

$$h_n^+\stackrel{~}{h}^+\mathrm{sinc}(\pi x_n)e^{i\pi x_n}\mathrm{where}x_n=f/f_mn,$$
(17)
and similarly for $`h_n^\times `$. Now suppose that the integration period $`T`$ is changed slightly to $`T^{}=T(1+ϵ)`$ where $`ϵ1`$. The Fourier coefficients of $`h(t)`$ evaluated over the interval $`T^{}`$ are related to those evaluated over the interval $`T`$ by

$$h_n^{}=\underset{j}{}h_j\mathrm{sinc}(\pi x_{jn})e^{i\pi x_{jn}},$$
(18)
where

$$x_{jn}=j\left(\frac{T^{}}{T}\right)n.$$
(19)
If the source frequency $`f`$ is an integer multiple of the new sample frequency $`\mathrm{\Delta }f^{}=1/T^{}`$, all the power will be concentrated in a single Fourier mode $`h_k^{}`$. Thus, the frequency spreading can be removed by finding the value of $`ϵ`$ that maximizes the spike in the stretched power spectra of $`h^+(t)`$ and $`h^\times (t)`$. This procedure yields an improved estimate of the source frequency, and a direct estimate of the Fourier amplitudes of the gravitational wave:

$$\stackrel{~}{h}^+h_{}^{+}{}_{k}{}^{}e^{i\delta }\mathrm{and}\stackrel{~}{h}^\times h_{}^{\times }{}_{k}{}^{}e^{i\delta }.$$
(20)
The quantity $`\delta =2\pi f(T^{}T)`$ accounts for the phase rotation imparted by the change in integration period. Once the Fourier amplitudes of the two polarizations have been estimated they can be used to calculate the amplitude $`𝒜`$, inclination $`\iota `$, polarization angle $`\psi `$ and orbital phase $`\phi _0`$ of the binary using the relations

$$\psi =\frac{1}{4}\left(\mathrm{arctan}\left(\frac{\stackrel{~}{h}_+^I+\stackrel{~}{h}_\times ^R}{\stackrel{~}{h}_\times ^I\stackrel{~}{h}_+^R}\right)+\mathrm{arctan}\left(\frac{\stackrel{~}{h}_+^I\stackrel{~}{h}_\times ^R}{\stackrel{~}{h}_+^R+\stackrel{~}{h}_\times ^I}\right)\right),$$
(21)
$$\phi _0=\frac{1}{2}\left(\mathrm{arctan}\left(\frac{\stackrel{~}{h}_+^I\stackrel{~}{h}_\times ^R}{\stackrel{~}{h}_+^R+\stackrel{~}{h}_\times ^I}\right)\mathrm{arctan}\left(\frac{\stackrel{~}{h}_+^I+\stackrel{~}{h}_\times ^R}{\stackrel{~}{h}_\times ^I\stackrel{~}{h}_+^R}\right)\right),$$
(22)
$$𝒜=\frac{A_++\sqrt{A_+^2A_\times ^2}}{2}$$
(23)
$$\iota =\mathrm{arccos}\left(\frac{A_\times }{A_++\sqrt{A_+^2A_\times ^2}}\right),$$
(24)
where

$`A_+={\displaystyle \frac{4(\mathrm{cos}2\psi \stackrel{~}{h}_+^R\mathrm{sin}2\psi \stackrel{~}{h}_\times ^R)}{\mathrm{cos}\phi _0}}`$
$`A_\times ={\displaystyle \frac{4(\mathrm{sin}2\psi \stackrel{~}{h}_+^R+\mathrm{cos}2\psi \stackrel{~}{h}_\times ^R)}{\mathrm{sin}\phi _0}},`$ (25)
and $`\stackrel{~}{h}_+^R=\mathrm{}(\stackrel{~}{h}^+)`$, $`\stackrel{~}{h}_+^I=\mathrm{}(\stackrel{~}{h}^+)`$, $`\stackrel{~}{h}_\times ^R=\mathrm{}(\stackrel{~}{h}^\times )`$, and $`\stackrel{~}{h}_\times ^I=\mathrm{}(\stackrel{~}{h}^\times )`$.

As a typical example of a LISA source, consider a White Dwarf binary with orbital period $`1000`$ seconds, comprised of two $`0.5M_{}`$ White Dwarfs, located 2 Kpc from Earth in the direction of the galactic center $`(\theta ,\varphi )=(95.54^{},266.83^{})`$. Setting $`\psi =111.2^{}`$, $`\iota =60^{}`$ and $`\phi _0=20.1^{}`$, and using the standard LISA noise curve, we arrive at the signal strengths shown in Fig. 1. The raw signal to noise $`\mathrm{SNR}_r`$, and the optimal signal to noise that can be achieved by matched filtering, $`\mathrm{SNR}_m`$, are defined:

$`\mathrm{SNR}_r=\left({\displaystyle \frac{_{\mathrm{bw}}S_h𝑑f}{_{\mathrm{bw}}S_n𝑑f}}\right)^{1/2}`$
$`\mathrm{SNR}_m=\left({\displaystyle _{\mathrm{bw}}}{\displaystyle \frac{2TS_hdf}{S_n}}\right)^{1/2},`$
where $`S_h`$ is the power spectral density of the signal and $`S_n`$ is the power spectral density of the noise, and the integration is across the bandwidth of the signal. Using these definitions we find that $`\mathrm{SNR}_r^\mathrm{I}=5.4`$, $`\mathrm{SNR}_r^{\mathrm{II}}=5.2`$, and the optimal signal to noise is $`\mathrm{SNR}_m=49.9`$. Applying the demodulation procedure yields the two gravitational wave polarizations shown in Fig. 2. The source parameters are recovered to an accuracy of $`\mathrm{\Delta }f=0.008f_m`$, $`\mathrm{\Delta }\theta =1.4^{}`$, $`\mathrm{\Delta }\varphi =0.2^{}`$, $`\mathrm{\Delta }𝒜/𝒜=0.054`$, $`\mathrm{\Delta }\iota =0.04^{}`$, $`\mathrm{\Delta }\psi =16.4^{}`$ and $`\mathrm{\Delta }\phi _0=8.6^{}`$. These uncertainties are comparable to those found using matched filtering, so the demodulation method is not only extremely fast, but also close to optimal.

The method continues to perform well when the SNR ratios are smaller, or when there are multiple sources in the same frequency band. Even with multiple sources the procedure delivers a good initial fit for the parameters describing each source. By combining the demodulation procedure with the gCLEAN algorithm described in Ref. cl2 , it will be possible to extract thousands of individual sources from the galactic background at a modest computational cost.

The simple power spike search will not work for rapidly evolving systems, such as supermassive black hole binaries, as even the demodulated power spectrum will be spread over thousands of frequency bins. If one thinks of a Fourier transform as a matched filter using sines and cosines, we see that we get a good match for monochromatic sources (the power spike), but a poor match for chirping sources. The solution is to employ a different kind of transform that uses filter functions that evolve in amplitude and frequency - the Fast Chirp Transform tom . Just as a Fast Fourier Transform (FFT) allows us to detect signals with a constant frequency, the Fast Chirp Transform (FCT) allows us to detect signals with variable frequency. By applying a FCT rather than a FFT to the demodulated signal, we can search for the sky locations that harbour the brightest chirping sources and extract their physical characteristics.

Much remains to be done in the field of LISA data analysis, but fast non-template based methods such as the total demodulation procedure and the Fast Chirp Transform will likely play a key role in future developments.

The work reported here is supported by the NASA EPSCoR program through Cooperative Agreement NCC5-579.


text:
# Neutrino flavour relaxation or neutrino oscillations?

Russian Research Centre ”Kurchatov Institute”,
123182, Kurchatov sq. 1, Moscow, Russia

We propose the new mechanism of neutrino flavour relaxation to explain the experimentally observed changes of initial neutrino flavour fluxes. The test of neutrino relaxation hypothesis is presented, using the data of modern reactor, solar and accelerator experiments. The final choice between the standard neutrino oscillations and the proposed neutrino flavour relaxation model can be done in future experiments.

PACS: 14.60.-z;14.60.Pq;76.20.+q

Now the phenomena of changes in initial neutrino flavour flux is observed in different neutrino experiments. The SNO experiment  evidently detects only 1/3 part of initial electron neutrino flux from the Sun and 2/3 part of muon and(or) tau neutrino fluxes. The Kamland reactor experiment  detected only 61% of expected electron antineutrino events from different reactors at mean distance of 180 km. Convincing evidence of initial neutrino flavour flux changes is observed also in Super-K  and MACRO  atmospheric neutrino data and K2K  accelerator experiment with muon neutrinos.

The standard common way of interpreting these results lies in neutrino oscillation hypothesis, first proposed by Bruno Pontecorvo  and developed in further works .

Here we discuss the alternative mechanism of neutrino flavour relaxation, which can also describe the observed changes of initial neutrino flavour fluxes with distance.

The proposed model is similar to the mechanism  of spin relaxation in random fluctuating magnetic field $`𝐁`$ with zero average $`<𝐁(t)>=0`$ and mean square fluctuating field value $`<𝐁^2(t)>0`$. The spin relaxation process is described by W.Pauli master equation .

Let us assume the existence of some small random fluctuating vacuum field $`\widehat{V}`$, causing the transitions between different lepton flavours. Such field $`\widehat{V}`$ can have mean zero value, with $`<\widehat{V}^2(t)>`$ different from 0 and $`(<\widehat{V}^2>)^{1/2}>m_{\nu _e},m_{\nu _\mu },m_{\nu _\tau }`$. Interactions of neutrino flavours with the vacuum field should lead to flavour relaxation process.

The time evolution of the neutrino states are governed by the Shrodinger equation:

$$i\frac{d}{dt}|\nu (t)>=\widehat{H}(t)|\nu (t)>,$$
(1)
where $`|\nu (t)>`$ is the neutrino vector of state and H(t) is the time-dependent Hamiltonian of the system, which form depends on in what basis it is given.

In flavour basis the total Hamiltonian in the random fluctuating vacuum field for rest reference frame can be written as

$$\widehat{H}_f(t)=\widehat{H}_m+\widehat{V}_f(t),$$
(2)
where

$$\widehat{H}_m=\left(\begin{array}{ccc}m_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right)\mathrm{and}\widehat{V}_f(t)=\left(\begin{array}{ccc}0& V_{e\mu }(t)& V_{e\tau }(t)\\ V_{\mu e}(t)& 0& V_{\mu \tau }(t)\\ V_{\tau e}(t)& V_{\tau \mu }(t)& 0\end{array}\right).$$
(3)
Here $`H_m`$ is the free Hamiltonian in mass basis. Note, that in this model the neutrino flavour states are assumed to be the same as their massive states (or in other words the mixing mass matrix is diagonal).

$`V_{\alpha \alpha ^{}}`$ $`(\alpha =e,\mu ,\tau )`$ is the vacuum field potential with mean value $`<V_{\alpha \alpha ^{}}(t)>=0`$ and $`<(V_{\alpha \alpha ^{}}(t))^2>0`$.

The neutrino flavour evolution in time can be obtained from Eq. (1-3), using the density matrix approach  and is given by Pauli master equation:

$$\frac{d}{dt}\left(\begin{array}{c}n_{\nu _e}\\ n_{\nu _\mu }\\ n_{\nu _\tau }\end{array}\right)=\left(\begin{array}{ccc}(W_{\mu e}+W_{\tau e})& W_{e\mu }& W_{e\tau }\\ W_{\mu e}& (W_{e\mu }+W_{\tau \mu })& W_{\mu \tau }\\ W_{\tau e}& W_{\tau \mu }& (W_{e\tau }+W_{\mu \tau })\end{array}\right)\left(\begin{array}{c}n_{\nu _e}\\ n_{\nu _\mu }\\ n_{\nu _\tau }\end{array}\right),$$
(4)
where $`n_\alpha (t)(\alpha =\nu _e,\nu _\mu ,\nu _\tau )`$ are the probability of observing neutrino with electron, muon or tau flavour, $`W_{\alpha \alpha \mathrm{`}}`$ corresponds to the neutrino transition rates from flavour $`\alpha \mathrm{`}`$ to flavour $`\alpha `$ and

$$\underset{\alpha =\nu _e,\nu _\mu ,\nu _\tau }{}n_\alpha (t)=1.$$

The general solution of Eq. (4) is given by the sum of two exponents and constant:

$$\left(\begin{array}{c}n_{\nu _e}(t)\\ n_{\nu _\mu }(t)\\ n_{\nu _\tau }(t)\end{array}\right)=\left(\begin{array}{c}a_1\\ a_2\\ a_3\end{array}\right)+\left(\begin{array}{c}b_1\\ b_2\\ b_3\end{array}\right)e^{t/T_1}+\left(\begin{array}{c}c_1\\ c_2\\ c_3\end{array}\right)e^{t/T_2},$$
(5)
with

$$\underset{i=1}{\overset{3}{}}a_i=1,\underset{i=1}{\overset{3}{}}b_i=0,\underset{i=1}{\overset{3}{}}c_i=0.$$
(6)
Note, that in relaxation model the neutrino flavour lepton numbers are violated, while the total sum of the numbers are constant (also as in the case of oscillation model).

To illustrate the possibilities of the relaxation model the particular ”simple” solution can be found under the following assumptions:

a) $`W_{e\mu }=W_{\mu e},W_{e\tau }=W_{\tau e},W_{\mu \tau }=W_{\tau \mu }`$,

b) $`W_{\mu \tau }W_{e\mu },W_{e\tau }`$

Now equation (5) looks as:

$$\left(\begin{array}{c}n_{\nu _e}(t)\\ n_{\nu _\mu }(t)\\ n_{\nu _\tau }(t)\end{array}\right)=\left(\begin{array}{c}1/3\\ 1/3\\ 1/3\end{array}\right)+b\left(\begin{array}{c}0\\ 1\\ 1\end{array}\right)e^{2W_{\mu \tau }t}+c\left(\begin{array}{c}1\\ 1/2\\ 1/2\end{array}\right)e^{3/2(W_{e\mu }+W_{e\tau })t}.$$
(7)
For small mass differences between different neutrino flavours ($`m_{\nu _e}m_{\nu _\mu }m_{\nu _\tau }`$ and $`\mathrm{\Delta }m_\nu m_\nu `$) with $`\gamma =E_\nu /m_\nu 1`$ the probability $`P`$ of observing the neutrino flavour with energy $`E_\nu `$ at distance $`r`$ from the source is given by:

$$\left(\begin{array}{c}P_{\nu _e}(r)\\ P_{\nu _\mu }(r)\\ P_{\nu _\tau }(r)\end{array}\right)=\left(\begin{array}{c}1/3\\ 1/3\\ 1/3\end{array}\right)+b\left(\begin{array}{c}0\\ 1\\ 1\end{array}\right)e^{r/(\mathrm{\Lambda }_0E_\nu )}+c\left(\begin{array}{c}1\\ 1/2\\ 1/2\end{array}\right)e^{r/(\mathrm{\Lambda }_1E_\nu )}.$$
(8)
where we use the new notation $`\mathrm{\Lambda }_0=(2W_{\mu \tau }m_\nu )^1`$

and $`\mathrm{\Lambda }_1=2/3[(W_{e\mu }+W_{e\tau })m_\nu ]^1`$.

Equation (8) can be solved for different initial experimental conditions of neutrino flavour fluxes.

For initial pure electron neutrino or antineutrino flux (like in the case of solar or reactor experiments) $`P_{\nu _e}(0)=1`$, $`P_{\nu _\mu }(0)=P_{\nu _\tau }(0)=0`$ the solution is given by:

$`P_{\nu _e}(r)`$ $`=`$ $`1/3[1+2e^{r/(\mathrm{\Lambda }_1E_\nu )}],`$
$`P_{\nu _\mu }(r)`$ $`=`$ $`1/3[1e^{r/(\mathrm{\Lambda }_1E_\nu )}],`$
$`P_{\nu _\tau }(r)`$ $`=`$ $`1/3[1e^{r/(\mathrm{\Lambda }_1E_\nu )}].`$ (9)
From Kamland experiment data  it is possible to estimate the value of $`\mathrm{\Lambda }_1`$ parameter. Taking the effective energy of reactor antineutrino flux to be equal to 4.8 MeV (here we take into account the threshold of Kamland detector $`E_{\stackrel{~}{\nu _e}}>`$ 3.4 MeV), mean reactor distance 180 km and $`R=P_{\stackrel{~}{\nu _e}(\mathrm{measured})}/P_{\stackrel{~}{\nu _e}(\mathrm{expected})}=0.61`$ we obtain

$$\mathrm{\Lambda }_1=43km/MeV.$$

Fig.1 illustrates the distance dependence of reactor antineutrino fluxes ( for $`E_{\stackrel{~}{\nu }}=4.8MeV`$).

Note, that the final flavour survival probability 1/3 is consistent with the SNO  experimental data with the ratio of measured neutrino flavour fluxes

$$\frac{CC}{NC}=0.306\pm 0.026(\mathrm{stat})\pm 0.024(\mathrm{syst}),$$
(10)
and with the absence of distortion in the measured neutrino spectrum at low energies, while such distortion is predicted by LMA solution  .

It agrees also with the Homestake results  for the ratio $`R`$ of observed and predicted by SSM  neutrino rates: $`R=0.34\pm 0.03`$.

The deviation from ”simple” relaxation model exists for Gallium experiments   with $`R=0.553\pm 0.034`$. But it is worth to mention, that the $`R`$ value depends on the accuracy of Standard Solar Model predictions for low energy fluxes.

For the case of pure initial muon neutrino or antineutrino flux (like in accelerator experiments) $`P_{\nu _\mu }(0)=1`$, $`P_{\nu _e}(0)=P_{\nu _\tau }(0)=0`$ the solution of Eq. (9) is:

$`P_{\nu _e}(r)`$ $`=`$ $`1/3[1e^{r/(\mathrm{\Lambda }_1E_\nu )}],`$
$`P_{\nu _\mu }(r)`$ $`=`$ $`1/3+1/2e^{r/(\mathrm{\Lambda }_0E_\nu )}+1/6e^{r/(\mathrm{\Lambda }_1E_\nu )},`$
$`P_{\nu _\tau }(r)`$ $`=`$ $`1/31/2e^{r/(\mathrm{\Lambda }_0E_\nu )}+1/6e^{r/(\mathrm{\Lambda }_1E_\nu )}.`$ (11)
Preliminary estimation of $`\mathrm{\Lambda }_0`$ parameter can be done using K2K accelerator experiment data . Taking the mean energy of muon neutrino flux $`E_\nu =1.3GeV`$, distance $`r=250km`$, $`R=P_{\nu _\mu (\mathrm{measured})}/P_{\nu _\mu (\mathrm{expected})}=0.70`$ we obtain:

$$\mathrm{\Lambda }_0=0.21km/MeV.$$

Fig.2 illustrates the distance dependence of neutrino fluxes from muon neutrino beam ($`E_\nu =1.3GeV`$).

It is necessary to note that the estimations of neutrino relaxation parameters $`\mathrm{\Lambda }_0`$ and $`\mathrm{\Lambda }_1`$ are very preliminary and were made mainly for illustration of relaxation model. More precise calculations can be done in future, using exact data of different experiments (taking also into account atmospheric neutrino data).

Finally we would like to emphasize the main difference between the standard oscillation theory and the proposed here relaxation model: the dependence of neutrino flavour fluxes from distance is described by the sum of constant and 2 exponents, instead of oscillation case. Fig.3 illustrates this difference for reactor antineutrino experiments.

Future reactor and accelerator neutrino experiments can provide the necessary data to choose between the neutrino oscillations and the proposed flavour relaxation model. Of course, the possibility to have a mixture of neutrino oscillation and relaxation models also exists.

Submitted to JETP Letters.