text:
# Two-component Bose-Einstein Condensates with Large Number of Vortices

## Abstract

We consider the condensate wavefunction of a rapidly rotating two-component Bose gas with an equal number of particles in each component. If the interactions between like and unlike species are very similar (as occurs for two hyperfine states of 87Rb or 23Na) we find that the two components contain identical rectangular vortex lattices, where the unit cell has an aspect ratio of $`\sqrt{3}`$, and one lattice is displaced to the center of the unit cell of the other. Our results are based on an exact evaluation of the vortex lattice energy in the large angular momentum (or quantum Hall) regime.

Experiments on rotating Bose gases have progressed rapidly in the last two years. Soon after the pioneer work at JILA  and ENS , the MIT group created a vortex lattice with as many as 160 vortices. Recently, the JILA group has invented an ingenious method to increase the angular momentum of a condensate by performing evaporative cooling on a rotating normal cloud. In this process, the system spins faster and faster as it is cooled, while remaining close to equilibrium. With such rapid progress, one expects that equilibrium Bose gases with even larger angular momenta may be produced in the near future.

At present, most experiments on vortex lattices are performed in single component Bose systems. It is natural to ask what happens in two-component Bose gases, such as those made up of two hyperfine spin states of the same atom. The vortex lattices in such systems are bound to be more intricate than those in single component condensates, as the vortices in different components can move relative to one another to minimize the energy. The purpose of this paper is to study the vortex lattices of two-component systems with large number of vortices, in what we call the “mean-field quantum Hall regime”. This is the regime where mean field theory remains valid so that each component (labeled by an index “$`i`$”, $`i=1,2`$) is characterized by a condensate wavefunction $`\mathrm{\Psi }_i`$; yet the angular momentum of the system is so high that $`\mathrm{\Psi }_i`$ is made up of the orbitals in the lowest Landau level in the plane perpendicular to the rotation axis. It has been shown recently that this regime will emerge in a three dimensional Bose gas at sufficiently high angular momenta . We focus on this regime because the wavefunction in this limit acquires an analytic structure which allows exact evaluation of the energy of a vortex lattice. As a result, it is possible to scan through a wide range of lattice structures which would be impractical for numerical calculations because of the time and the accuracy required. Although not directly applicable to current experiments on vortex lattices (which are performed at lower angular momenta), the physics of the mean field quantum Hall regime is still quite relevant as the vortex lattices in these two regimes are connected continuously to each other.

One special feature of the majority of two-component gases so far studied, (notably mixtures of hyperfine states of 87Rb  in magnetic traps, or 23Na  in optical traps) is that the interactions between like species (denoted $`g_1`$ and $`g_2`$) and unlike species (denoted $`g_{12}`$) are very similar, within a few percent of each other. Thus, if there are an equal number of bosons in each component, and if the trapping potentials of the two components are made identical, then the two components will be the same size and contain the same density of vortices. In this case, one expects that each component will contain identical vortex lattices, with one lattice displaced relative to the other. While we are mainly interested in the experimentally relevant cases, where $`g_1g_2g_{12}`$, considerable insight is gained by studying vortex lattices as a function of the interactions. Considering the case $`g_1g_2g_{12}`$, we find a wide range of vortex lattice structures as the parameter $`\alpha =g_{12}/\sqrt{g_1g_2}`$ is varied. The vortex lattice has a fixed structure over certain intervals of $`\alpha `$, while it varies continuously in others. The structure near the isotropic point $`g_1=g_2=g_{12}`$ consists of identical rectangular lattices in both components, with one displaced to the center of the unit cell of the other. The aspect ratio of the unit cell changes with $`\alpha `$, and is $`\sqrt{3}`$ when $`\alpha =1`$.

The mean-field quantum Hall regime: The condensate wavefunctions $`\mathrm{\Psi }_1`$ and $`\mathrm{\Psi }_2`$ of a two-component rotating Bose gas are determined by minimizing the grand potential $`K=E\mathrm{\Omega }L_z\mu _1N_1\mu _2N_2`$, where $`E`$ is the energy of the system, $`\mathrm{\Omega }`$ is the rotational frequency, $`L_z`$ is the angular momentum along $`z`$, and $`\mu _i`$ ($`i=1,2`$) are the chemical potentials fixing the number of bosons $`N_1`$ and $`N_2`$ in each component. For simplicity, we assume identical trapping potentials for each component. We consider a cigar shaped trap with the symmetry axis $`z`$ coinciding with the axis of rotation. As discuss in , the slow variation of the trapping potential along $`z`$ allows one to apply a Thomas-Fermi approximation for the $`z`$ dependence of $`\mathrm{\Psi }_i`$ and write $`K`$ as $`dzd𝐫𝒦(𝐫,z)`$,

$`𝒦(𝐫,z)={\displaystyle \underset{i=1,2}{}}\mathrm{\Psi }_i^{}[h\mu _i(z)]\mathrm{\Psi }_i+𝒱,`$ (1)
$`h={\displaystyle \frac{1}{2M}}\left({\displaystyle \frac{\mathrm{}}{i}}M\mathrm{\Omega }\widehat{𝐳}\times 𝐫\right)^2+{\displaystyle \frac{1}{2}}M(\omega _{}^2\mathrm{\Omega }^2)r^2,`$ (2)
$`𝒱={\displaystyle \frac{1}{2}}g_1|\mathrm{\Psi }_1|^4+{\displaystyle \frac{1}{2}}g_2|\mathrm{\Psi }_2|^4+g_{12}|\mathrm{\Psi }_1|^2|\mathrm{\Psi }_2|^2,`$ (3)
with $`𝐫=(x,y)`$, $`\mu _i(z)=\mu _i\frac{1}{2}M\omega _z^2z^2`$, $`g_i=4\pi \mathrm{}^2a_i/M`$, $`(i=1,2)`$, and $`g_{12}=4\pi \mathrm{}^2a_{12}/M`$, where $`a_i`$, and $`a_{ij}`$ are the s-wave scattering lengths between like and unlike bosons respectively. As $`z`$ is treated as a parameter, it is convenient to write $`\mathrm{\Psi }_i=\sqrt{n_i(z)}\mathrm{\Phi }_i(𝐫;z)`$, with $`|\mathrm{\Phi }_i(𝐫;z)|^2d^2r=1`$. The number constraint $`d𝐫dz|\mathrm{\Psi }|^2=N_i`$ becomes $`n_i(z)dz=N_i`$.

As $`\mathrm{\Omega }`$ approaches $`\omega _{}`$, the wavefunctions $`\mathrm{\Phi }_i`$ are made up of the orbitals $`u_m(𝐫)`$ in the lowest Landau level in the $`xy`$-plane, $`\mathrm{\Phi }_i(𝐫,z)=_{m=0}^{\mathrm{}}c_m(z)u_m(𝐫)`$, where $`u_m(𝐫)=(2\pi m!)^{1/2}[(x+iy)/d]^me^{r^2/2d^2}`$, and $`d=\sqrt{\mathrm{}/M\omega _{}}`$. The potential $`𝒦`$ then becomes

$$𝒦=\underset{i=1,2}{}\left[\mathrm{}(\omega _{}\mathrm{\Omega })\frac{r^2_i}{d^2}\mu _i(z)+\mathrm{}\omega _{}\right]n_i(z)+𝒱$$
(4)
where $`r^2_i=r^2|\mathrm{\Phi }_i|^2d^2r`$, and

$$𝒱=\mathrm{d}^2r\left[\frac{1}{2}\underset{i=1,2}{}g_in_i^2|\mathrm{\Phi }_i|^4+g_{12}n_1n_2|\mathrm{\Phi }_1|^2|\mathrm{\Phi }_2|^2\right].$$
(5)
As shown in ref., wavefunctions in the lowest Landau level (not normalized) can be written as

$$\varphi (𝐫)=\lambda \underset{\alpha }{}(wa_\alpha )e^{r^2/2d^2},w=x+iy.$$
(6)
where $`\lambda `$ is an arbitrary constant and $`\{a_\alpha \}`$ are the zeros of $`\varphi `$. If the zeros form a infinite lattice with unit cell size $`v_c`$, it is shown in  that $`|\varphi |^2`$ is a product of a Gaussian and a function periodic under lattice translation. i.e.

$$|\varphi |^2=e^{r^2/\sigma ^2}g(𝐫),g(𝐫)=g(𝐫+𝐑),$$
(7)
where $`𝐑=n_1𝐁_1+n_2𝐁_2`$, $`n_i`$ are integers, and $`𝐁_1`$, $`𝐁_2`$ are the basis vectors of the lattice. The width $`\sigma `$ reflects the number of vortices of the system. It is given by

$$\sigma ^2=d^2\pi v_c^1.$$
(8)
The periodicity of $`g(𝐫)`$ implies $`g(𝐫)=v_c^1_𝐊g_𝐊e^{i𝐊𝐫}`$, where $`\{𝐊\}`$ are the reciprocal lattice vectors.

In the following, we shall consider a two-component Bose gas with equal particle numbers and trapping potentials, and with interactions $`g_1g_2g_{12}`$. If $`g_1=g_2`$, the two components are identical and one expects that each will contain identical vortex lattices, translated with respect to one another. Sufficiently small differences in $`g_1g_2`$ should not change this structure (though changes may occur in the density profiles $`n_i(z)`$, the parameters of the lattice, and the relative displacement $`r_0`$). This structure persists because, even when $`g_1g_2`$, the two components contain equal vorticity, hence an equal density of vortices. The potential energy is minimized by interlacing the two lattices; if the vortex lattice in one component were to deform, the other has to follow to keep the interaction energy at a minimum. We shall therefore consider normalized condensates $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ in Eq. (4) with densities of the form

$`|\mathrm{\Phi }_1|^2`$ $`=`$ $`(\pi \sigma ^2)^1_𝐊\stackrel{~}{g}_𝐊e^{i𝐊𝐫}e^{r^2/\sigma ^2}`$ (9)
$`|\mathrm{\Phi }_2|^2`$ $`=`$ $`(\pi \sigma ^2)^1_𝐊\stackrel{~}{g}_𝐊e^{i𝐊(𝐫𝐫_o)}e^{r^2/\sigma ^2}`$ (10)
$`\stackrel{~}{g}_𝐊`$ $`=`$ $`g_𝐊/\left(_𝐊^{}g_𝐊^{}e^{\sigma ^2𝐊_{}^{}{}_{}{}^{2}/4}\right).`$ (11)
The wavefunction is described by variational parameters $`n_i(z)`$, $`\sigma ^2`$, the basis vectors $`𝐁_i`$ (which determine the unit cell size $`v_c`$), and the relative displacement $`𝐫_o`$.

By integrating Eqs. (9) and (10), one sees that up to terms of relative order $`v_c/\sigma ^2`$ the cloud’s radius is $`r^2_1=r^2_2=\sigma ^2`$. Defining the quantities $`I`$ and $`I_{12}`$ as $`|\mathrm{\Phi }_i|^4\mathrm{d}^2rI/(\pi \sigma ^2)`$ and $`|\mathrm{\Phi }_1|^2|\mathrm{\Phi }_2|^2\mathrm{d}^2r=I_{12}/(\pi \sigma ^2)`$, we have

$`I`$ $`=`$ $`{\displaystyle \underset{𝐊,𝐊^{}}{}}\stackrel{~}{g}_𝐊\stackrel{~}{g}_𝐊^{}e^{\sigma ^2|𝐊+𝐊^{}|^2/4},`$ (12)
$`I_{12}`$ $`=`$ $`{\displaystyle \underset{𝐊}{}}\stackrel{~}{g}_𝐊\stackrel{~}{g}_𝐊^{}e^{i𝐊^{}𝐫_o}e^{\sigma ^2|𝐊+𝐊^{}|^2/4},`$ (13)
and the potential $`𝒦`$ takes the form

$`𝒦=`$ $`\left[\mu (z)\mathrm{}\omega _{}\mathrm{}(\omega _{}\mathrm{\Omega })(\sigma ^2/d^2)\right](n_1+n_2)`$ (14)
$`+(d^2/2\sigma ^2)\left(n_1^2g_1I+n_2^2g_2I+2n_1n_2g_{12}I_{12}\right).`$ (15)
To evaluate $`𝒦`$, we need an explicit expression for the coefficients $`g_𝐊`$. These coefficients are most easily derived by introducing the complex representation for the basis vectors, $`b_i(\widehat{𝐱}+i\widehat{𝐲})𝐁_i`$. The area of the unit cell is then $`v_c=i(b_1^{}b_2b_2^{}b_1)/2`$. If we orient the lattice so that $`b_1`$ is real, i.e. $`𝐁_1=b_1\widehat{𝐱}`$, $`𝐁_2=b_1(u\widehat{𝐱}+v\widehat{𝐲})`$, we then have

$$b_2b_1(u+iv),v_c=b_1^2v.$$
(16)
In the Appendix, we show that a function $`\mathrm{\Phi }`$ in the lowest Landau level describing a regular vortex lattice contained in a cylindrically symmetric cloud will have the form $`\varphi (𝐫)=f(w)e^{r^2/2d^2}`$, with $`wx+iy`$, and

$$f(w)=\theta (\zeta ,\tau )e^{\pi w^2/2v_c}$$
(17)
where $`\zeta =w/b_1=(x+iy)/b_1`$, $`\tau =u+iv=b_2/b_1`$, and $`\theta `$ is the Jacobi theta-function defined as

$$\theta (\zeta ,\tau )=\frac{1}{i}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}(1)^ne^{i\pi \tau (n+1/2)^2}e^{2\pi i\zeta (n+1/2)}.$$
(18)
The density $`|\varphi |^2`$ is therefore of the form Eq. (7), with $`\sigma `$ given in Eq. (8), and

$$g(𝐫)=|\theta (\zeta ,\tau )\mathrm{e}\mathrm{x}\mathrm{p}(\pi y^2/v_c)|^2.$$
(19)
The Jacobi theta-function has the quasi-periodic properties

$`\theta (\zeta +1,\tau )`$ $`=`$ $`\theta (\zeta ,\tau ),`$ (20)
$`\theta (\zeta +\tau ,\tau )`$ $`=`$ $`e^{i\pi (\tau +2\zeta )}\theta (\zeta ,\tau ),`$ (21)
which implies the periodic property $`g(𝐫)=g(𝐫+𝐑)`$. The Fourier coefficients of $`g(𝐫)`$ are

$$g_𝐊=(1)^{m_1+m_2+m_1m_2}e^{v_c|𝐊|^2/8\pi }\sqrt{\frac{v_c}{2}},$$
(22)
where $`𝐊=m_1𝐊_1+m_2𝐊_2`$, and $`𝐊_i`$ are the basis vector of the reciprocal lattice, $`𝐊_1=(2\pi /v_c)𝐁_2\times \widehat{𝐳}`$, $`𝐊_2=(2\pi /v_c)\widehat{𝐳}\times 𝐁_1`$, and

$$v_c𝐊^2=\frac{(2\pi )^2}{v}\left((vm_1)^2+(m_2um_1)^2\right)$$
(23)
Since we work in the limit of large vortex number, the size of the cloud is much larger than the unit cell, i.e. $`\pi \sigma ^2/v_c>>1`$. We can therefore ignore all $`𝐊+𝐊^{}0`$ terms in Eqs. (12) and (13), since $`\sigma ^2𝐊^2>\pi \sigma ^2/v_c`$. We then have

$$I=\underset{𝐊}{}\left|\frac{g_𝐊}{g_\mathrm{𝟎}}\right|^2,I_{12}=\underset{𝐊}{}\left|\frac{g_𝐊}{g_\mathrm{𝟎}}\right|^2\mathrm{cos}𝐊𝐫_o,$$
(24)
where $`g_𝐊`$ is given by Eq. (22) and $`𝐊`$-sum is over the integers $`m_1,m_2`$. Since the expressions of $`I`$ and $`I_{12}`$ in Eq. (24) are independent of $`\sigma ^2`$, the minimization of $`𝒦`$ in Eq. (15) with respect to $`\sigma ^2`$ and $`n_i`$ become very simple. The optimum $`\sigma ^2,𝐧=(n_1,n_2)`$, and $`𝒦`$ are given by

$`\sigma ^2`$ $`=`$ $`d^2\left(\mu (z)\mathrm{}\omega _{}\right)/\left(3\mathrm{}(\omega _{}\mathrm{\Omega })\right),`$ (25)
$`𝐧(z)`$ $`=`$ $`(2/3)(\sigma ^2/d^2)(\mu (z)\mathrm{}\omega _{})𝐆^1\mathrm{𝟏},`$ (26)
$`𝒦`$ $`=`$ $`(1/3)(\mu (z)\mathrm{}\omega _{})\mathbf{\hspace{0.17em}\hspace{0.17em}1}𝐧(z),`$ (27)
$`𝐆`$ $`=`$ $`\left(\begin{array}{cc}g_1I& g_{12}I_{12}\\ g_{12}I_{12}& g_2I\end{array}\right),\mathbf{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}=\left(\begin{array}{c}1\\ 1\end{array}\right),`$ (32)
It is clear from Eqs. (25) through (32) that the solution for the case where $`g_1g_2<<|g_1+g_2|`$ is very close to that of $`g_1=g_2`$. The lattice shape (parameterized by $`𝐫_0`$, $`u`$, and $`v`$) enters the grand potential only through the factor $`\mathrm{𝟏}𝐆^\mathrm{𝟏}\mathrm{𝟏}`$. When $`g_1=g_2`$ this term is inversely proportional to $`J=I+\alpha I_{12}`$, and the most favorable lattice is the one that minimizes $`J`$.

Summary of Results: It is interesting to compare the two-component case with the single-component case. In the latter system, energy minimization reduces to minimizing $`I`$. The only local minimum is the triangular lattice, where $`I=1.1596`$; the square lattice is a saddle point with $`I=1.1803`$. The minute difference between these values of $`I`$ makes a simple numerical minimization of (1) challenging and illustrates the utility of the analytic scheme used here.

For a two-component Bose gas, the most favorable lattice minimizes $`I+\alpha I_{12}`$. In the minimization it is convenient to measure lengths in unit of the the basis vector $`𝐁_1=b_1\widehat{𝐱}`$, and write complex representation of $`𝐁_2`$ and $`𝐫_o`$ as $`\tau =u+iv=|\tau |e^{i\eta }`$ and $`r_oa+b\tau `$ respectively. The phase diagram of the vortex lattice as a function of the ratio $`\alpha =g_{12}/g`$ is shown in Figs. 1 and 2. The major features are

(a) $`\alpha <0`$: In this region the vortices of the two components coincide with each other ($`a=b=0`$) to form a triangular lattice ($`\tau =e^{i\pi /3}`$).

(b) $`0<\alpha <0.172`$: the vortex lattice in each component remains triangular. However, $`r_o`$ undergoes a first order change so that one lattice is displaced to the center of the triangle of the other ($`a=b=1/3`$). The lattice type (characterized by $`\tau =e^{i\pi /3}`$) remains constant within this interval.

(c) $`0.172<\alpha =0.373`$: At $`\alpha =0.172`$, $`𝐫_o`$ jumps from the center of the triangle (i.e. half of the unit cell) to the center of the (rhombic) unit cell ($`a=b=1/2`$). The angle $`\eta `$ jumps from $`60^o`$ to $`67.95^o`$ at $`\alpha =0.172`$, and increases continuously to $`90^o`$ as $`\alpha `$ increases to $`0.372`$. As a result, the lattice type is no longer fixed and the unit cell is a rhombus. The modulus of $`\tau `$, however, remains fixed across this region.

(d) $`0.373<\alpha <0.926`$: The two lattices are “mode-locked” into a centered square structure throughout the entire interval $`(\tau =i,a=b=1/2)`$.

(e) $`0.926<\alpha `$: The lattice type again varies continuously with interaction $`\alpha `$. Each component’s vortex lattice has a rectangular unit cell ($`\eta =\pi /2`$) whose aspect ratio $`|\tau |`$ increases with $`\alpha `$. Both 87Rb and 23Na have interaction parameters within this range. At $`\alpha =1`$, ($`g_1=g_2=g_{12}=g`$), the aspect ratio is $`\sqrt{3}`$. If one ignores the difference between the components, the combined lattice is triangular, as is expected.

It is interesting to note that in the absence of rotation, the two components change from miscible to immiscible when $`\alpha `$ increase beyond $`1`$. No such change, however, happens at $`\alpha =1`$ in the high angular momentum limit. This qualitative difference in behavior occurs because the presence of a vortex lattice naturally modulates the density of each component, with the high density regions of one fluid coincident with the low density regions of the other. Thus the system is effectively phase separated whenever staggered vortex lattices are present, even for $`\alpha <1`$. In particular, the vortex lattice near $`\alpha =1`$ (above or below) is made up of alternating rows of vortices of each component, (see Fig. 1), and the system therefore contains stripes in which one component has a high density and the other component has a very low density. As $`\alpha `$ increases, the stripes become more pronounced.

Final Remarks: The diversity of the vortex lattice structures in the two-component Bose gas has once again demonstrated the rich properties of these systems. Our calculation, based on exact evaluation of the vortex energy, assumes a perfect lattice. Considering the long relaxation times in clouds of dilute atoms, one might see more complicated structures, where patches of vortex domains are separated by defects or grain boundaries. Nevertheless, the underlying equilibrium structure should be reflected within each vortex domain.

So far, we have only discussed two-components systems with simple interpenetrating Bravais lattices. Our method is more general in that it allows the exact evaluation of the energy of an arbitrary regular lattice (with however complicated unit cell decoration). Such structures may be favored when the particles in each component have different numbers, trapping potentials, or masses (as in the case of 23Na-87Rb mixtures).

Appendix:  The general form of a vortex lattice in the lowest Landau level is $`\mathrm{\Psi }(x,y)=f(w)e^{r^2/2d^2}`$, where $`w=x+iy`$ and $`f`$ is an entire function whose zeros form a regular lattice $`\{b=nb_1+n_2b_2\}`$, where $`n_i`$ are integers, and $`b_1`$ and $`b_2`$ are the complex basis vectors, ($`b_2=b_1(u+iv)`$). Since the Jacobi theta function $`\theta ((x+iy)/b_1,u+iv)`$ is an entire function with exactly these zeros, we have $`f(w)=\theta (\zeta ,\tau )h(\zeta )`$, where $`\zeta =(x+iy)/b_1=\overline{x}+i\overline{y}`$, $`\tau =u+iv=b_2/b_1`$, and $`h(\zeta )`$ is an entire function without zeros. To ensure the normalizability of $`\mathrm{\Psi }`$, this function can only be of the form $`h(\zeta )=\mathrm{exp}(c_1\zeta +c_2\zeta ^2)`$. It is straightforward to show that

$`|\theta (\zeta ,\tau )|^2={\displaystyle \underset{m}{}}(1)^me^{2\pi im\overline{x}}e^{\pi vm^2/2}L_m`$ (33)
$`L_m={\displaystyle \frac{1}{2}}{\displaystyle \underset{m^{}}{}}\left(1e^{i\pi (m+m^{})}\right)e^{\left(i\pi um2\pi \overline{y}\pi vm^{}/2\right)m^{}}`$ (34)
By applying the Poisson summation formula to to Eq. (34), we have

$$L_m=\sqrt{\frac{2}{v}}\underset{k}{}(1)^{(m+1)k}e^{\left(\pi (k+um+2i\overline{y})^2/2v\right)}.$$
(35)
We thus have

$$|\theta (\zeta ,\tau )|^2=\left[\frac{1}{v_c}\underset{𝐊}{}g_𝐊e^{i𝐊𝐫}\right]e^{2\pi y^2/v_c}$$
(36)
where $`𝐫=x\widehat{x}+y\widehat{y}`$, $`𝐊=(2\pi m\widehat{x}2\pi (n+um)/v\widehat{y})/b_1`$, and $`g_𝐊`$ is given by (22). The density of the system is then $`|\mathrm{\Psi }(𝐫)|^2=|\theta (\zeta ,\tau )|^2|e^{c_1\zeta +c_2\zeta ^2}|^2e^{r^2/d^2}`$. For a vortex lattice with inversion symmetry about the origin $`𝐫=0`$, we have $`c_1=0`$. In addition, if the cloud’s envelope is cylindrically symmetric, we have $`c_2=\pi /(2v_c)`$, which gives Eqs. (7), (8), (17) and (19). Similar approaches have been used by Tkachenko  and Abrikosov  in their respective studies of 4He and type two superconductors.

This work is supported by NASA Grants NAG8-1441, NAG8-1765, and by NSF Grants DMR-0109255, DMR-0071630.


format: text

text:
# 1 Introduction

## 1 Introduction

In geophysics the matter density profile of the Earth is determined by measurements of the seismic wave velocity profile. Seismic method is used widely also for the oil search, in particular, in the decision process before a trial well is drilled. This method has, however, shortcomings and alternatives are highly appreciated.

Enormous penetration ability of neutrinos is a source of strong temptation to use them in tomography and geological studies ,,, . The first proposal was based on inelastic scattering of high energy neutrinos produced by accelerators . The authors have suggested several methods which use difference of the neutrino interactions in matter with different chemical composition and density.

Qualitatively different proposal is based on the elastic forward scattering (refraction) of neutrinos ,,. Refraction modifies properties of neutrino oscillations changing both the oscillation length and depth . The effect depends on density of medium. Consequently, an appearance of layers with different densities on the way of neutrinos changes the oscillation pattern. In this connection, some possibilities to use the neutrino superbeams as well as the beams from neutrino factories have been considered .

There is a number of studies of the Earth matter effect on oscillations of solar neutrinos. It was marked that the solar neutrinos can give information about large scale density distribution inside the Earth. In particular, properties of the mantle and the core of the Earth can be studied .

The Earth matter effect has been considered for the $`{}_{}{}^{7}Be`$ neutrinos  in connection to BOREXINO  and KamLAND  experiments. It was found that the effect is negligible for the LMA solution, at least for these detectors.

Main conditions to use neutrinos for the geological studies and searches for oil and minerals are:

1) sensitivity to small scale structures: $`d`$ 100 km;

2) possibility to move both a source and detector of neutrinos, so that substantial part of outer layers of the Earth mantle can be scanned.

No realistic proposals which satisfy these conditions have been published so far. In this paper we will consider a possibility which is based on:

* Detection of the solar $`{}_{}{}^{7}Be`$ neutrinos: For the LMA solution of the solar neutrino problem the oscillation length of neutrinos with $`E_{Be}`$ 0.86 MeV is in the required range lν/2 $``$ (10 $`÷`$ 20) km.
* Measuring the time dependence of the neutrino flux by deep underwater detector - submarine.

The paper is organized as follows. In section 2 we consider the effect of thin layers of matter on neutrino oscillations (see  for some related work). We study dependence of the effect on neutrino parameters and properties of layers. In section 3 we consider a possibility to identify the effect of cavities. In section 4 we discuss problems of detections. Concluding remarks are given in section 5.

## 2 Effects of thin layers on neutrino oscillations

Let us consider a system of two mixed neutrinos with values of oscillation parameters $`\mathrm{\Delta }m^2`$ and $`\mathrm{sin}^22\theta `$ from the LMA solution region . (We will comment on the effect of the third neutrino at the end of this section). In the case of the LMA solution, the electron neutrino produced in the center of the sun is adiabatically converted to a combination of the mass eigenstates $`\nu _1`$ and $`\nu _2`$ which is determined by the mixing angle, $`\theta _m^0`$, in the production point:

$$\nu _e\mathrm{cos}\theta _m^0\nu _1+\mathrm{sin}\theta _m^0\nu _2.$$
(1)
The angle $`\theta _m^0`$ is given by:

$$\mathrm{tan}^22\theta _m^0=\mathrm{tan}^22\theta \left(1\frac{2V_e^0E}{\mathrm{\Delta }m^2\mathrm{cos}2\theta }\right)^2,$$
(2)
where $`E`$ is the neutrino energy, $`V_e^0`$ is the matter potential for neutrinos in the production point. $`V_e=\sqrt{2}G_F\rho N_AY_e`$, $`G_F`$ is the Fermi coupling constant, $`\rho `$ is the density of medium, $`N_A`$ is the Avogadro number and $`Y_e`$ is the number of electrons per nucleon; $`V_e^0=V_e(\rho ^0,Y_e^0)`$.

Conversion effect should be averaged over the neutrino production region, and in what follows we will describe this averaging by the effective mixing angle $`\overline{\theta }_m^0`$.

Due to loss of coherence , neutrinos arrive at the surface of the Earth as incoherent fluxes of $`\nu _1`$ and $`\nu _2`$ with relative admixtures given by $`\mathrm{cos}^2\overline{\theta }_m^0`$ and $`\mathrm{sin}^2\overline{\theta }_m^0`$ correspondingly.

Let us consider oscillations inside the Earth. The probability to find $`\nu _e`$ in the detector can be written as:

$$P=\mathrm{cos}^2\overline{\theta }_m^0P_{1e}+\mathrm{sin}^2\overline{\theta }_m^0P_{2e}=\mathrm{cos}2\overline{\theta }_m^0P_{1e}+\mathrm{sin}^2\overline{\theta }_m^0,$$
(3)
where $`P_{1e}`$ and $`P_{2e}`$ are the probabilities of $`\nu _1\nu _e`$ and $`\nu _2\nu _e`$ transitions in matter of the Earth correspondingly.

Suppose neutrinos cross the cavity with length $`d`$ and density $`\rho _d`$ which differs from the average density of surrounding matter, $`\rho `$. That is, the neutrinos cross consecutively three layers: with lengths and densities ($`l_1`$, $`\rho `$), ($`d`$, $`\rho _d`$), ($`l_2`$, $`\rho `$). The total length of the trajectory is $`L=l_1+d+l_2`$.

Let us find an effect of the cavity on oscillations. Introducing, $`P_{1e}^0`$, the $`\nu _1\nu _e`$ oscillation probability in absence of cavity (it corresponds to oscillations in the unique layer of density $`\rho `$ and length $`L`$) we can write $`P`$ as:

$$P=P^0+\mathrm{\Delta }P,$$
(4)
$$P^0=\mathrm{sin}^2\overline{\theta }_m^0+\mathrm{cos}2\overline{\theta }_m^0P_{1e}^0,$$
(5)
$$\mathrm{\Delta }P=\mathrm{cos}2\overline{\theta }_m^0(P_{1e}P_{1e}^0).$$
(6)
Here $`P^0`$ is the probability to find $`\nu _e`$ in the absence of cavity and $`\mathrm{\Delta }P`$ is the change of the probability due to the cavity effect.

We will calculate the effect for the beryllium neutrino line with energy $`E_{Be}=0.86`$ MeV. For such an energy the matter effect on mixing and oscillations is very weak. It is determined by a small parameter:

$$ϵ\frac{2V_eE}{\mathrm{\Delta }m^2}2.610^3\left(\frac{\rho }{2.7\mathrm{g}/\mathrm{cm}^3}\right)\left(\frac{710^5\mathrm{eV}^2}{\mathrm{\Delta }m^2}\right)\left(\frac{Y_e}{0.5}\right)$$
(7)
which characterizes deviations of the mixing angle and the oscillation length in medium from their vacuum values. We will use $`ϵ`$ as an expansion parameter and find effects in the lowest order in $`ϵ`$. The oscillation length is determined by $`\mathrm{\Delta }m^2`$:

$$l_\nu ^ml_\nu 30.7\mathrm{km}\left(\frac{710^5\mathrm{eV}^2}{\mathrm{\Delta }m^2}\right).$$
(8)
The probability of $`\nu _1\nu _e`$ oscillations equals

$$P_{1e}=|\mathrm{cos}\theta S_{11}+\mathrm{sin}\theta S_{12}|^2,$$
(9)
where $`S=S_{ij}`$ is the evolution matrix in the mass eigenstate basis, $`\stackrel{}{\nu }`$. The matrix $`S`$ can be written as:

$$S=UD_2U^{}U_dD_dU_d^{}UD_1U^{}.$$
(10)
Here $`UU(\rho )`$ and $`U_dU(\rho _d)`$ are the mixing matrices of mass eigenstates in the layers 1, 2 and in the cavity correspondingly. They are defined in such a way that $`\stackrel{}{\nu }=U\stackrel{}{\nu }_m`$, where $`\stackrel{}{\nu }_m`$ is the basis of eigenstates of the Hamiltonean in matter. In (10), the diagonal matrices,

$$D_idiag(1,e^{i\mathrm{\Phi }_i}),i=1,2,d,$$
(11)
describe evolution of the neutrino eigenstates, $`\stackrel{}{\nu }_m`$, in layer $`i`$. $`\mathrm{\Phi }_i`$ is the phase of oscillations:

$$\mathrm{\Phi }_i=\frac{\mathrm{\Delta }m^2l_i}{2E}\sqrt{(\mathrm{cos}2\theta ϵ_i)^2+\mathrm{sin}^22\theta },$$
(12)
$`ϵ_iϵ(\rho _i)`$. In absence of cavity we have

$$S=S_0=UDU^+,$$
(13)
where $`D=diag(1,e^{i\mathrm{\Phi }})`$ and

$$\mathrm{\Phi }=\frac{\mathrm{\Delta }m^2L}{2E}\sqrt{(\mathrm{cos}2\theta ϵ)^2+\mathrm{sin}^22\theta }$$
(14)
is the total oscillation phase in absence of the cavity.

Using Eqs. (9,10,11), the straightforward calculation leads to the $`\nu _1\nu _e`$ probability in the lowest order in $`ϵ`$:

$$P_{1e}P_{1e}^0+ϵ\xi \mathrm{sin}^22\theta \mathrm{sin}\left(\mathrm{\Phi }_2+\frac{\mathrm{\Phi }_d}{2}\right)\mathrm{sin}\frac{\mathrm{\Phi }_d}{2},$$
(15)
where

$$\xi \frac{(Y_e\rho )(Y_{ed}\rho _d)}{Y_e\rho }$$
(16)
is the density contrast. In absence of cavity the probability equals

$$P_{1e}^0=\mathrm{cos}^2\theta ϵ\mathrm{sin}^22\theta \mathrm{sin}^2\frac{\mathrm{\Phi }}{2}.$$
(17)
And according to eq. (6) and (15) the effect of cavity can be written as

$$\mathrm{\Delta }P=ϵ\xi \mathrm{cos}2\overline{\theta }_m^0\mathrm{sin}^22\theta \mathrm{sin}\left(\mathrm{\Phi }_2+\frac{\mathrm{\Phi }_d}{2}\right)\mathrm{sin}\frac{\mathrm{\Phi }_d}{2}.$$
(18)
The following remarks are in order:

According to eq. (18) the effect of cavity is proportional to expansion parameter $`ϵ`$, the density contrast, $`\xi `$, and $`\mathrm{cos}2\overline{\theta }_m^0`$. For maximal mixing in the production region the effect would be zero. For the best fit point of the LMA solution we get $`\mathrm{cos}2\overline{\theta }_m^00.3`$.

Important feature of the result (18) is that $`\mathrm{\Delta }P`$ does not depend on the oscillation effect in the first layer (before cavity), in particular, it does not depend on the phase $`\mathrm{\Phi }_1`$. In contrast, $`\mathrm{\Delta }P`$ does depend on $`\mathrm{\Phi }_2`$ – the oscillation phase in the layer between the cavity and detector. This property appears in the lowest order in $`ϵ`$ and is related to the fact that initial state is the mass eigenstate and final state (at the detector) is the flavor state. It can be shown that if initial state is $`\nu _e`$ and final state – $`\nu _1`$, the probability $`\mathrm{\Delta }P`$ does not depend on $`\mathrm{\Phi }_2`$ (phase between cavity and detector) but will depend on $`\mathrm{\Phi }_1`$ the phase acquired in the first layer. Thus, the interchange of $`\nu _1`$ and $`\nu _e`$ in the initial and final states leads to change $`\mathrm{\Phi }_2\mathrm{\Phi }_1`$. This feature is important for identification of the cavity effect (see section 3).

Let us estimate the size of possible effect depending on parameters of neutrinos and the cavity. The cavity produces the change of the event rate during the time when the neutrino beam crosses the cavity. The relative value of the change $`\mathrm{\Delta }N/N=\mathrm{\Delta }P/P`$ equals in the lowest order in $`ϵ`$

$$\frac{\mathrm{\Delta }P}{P^0}=f(\mathrm{\Delta }m^2,\theta )ϵ\xi \mathrm{sin}\left(\mathrm{\Phi }_2+\frac{\mathrm{\Phi }_d}{2}\right)\mathrm{sin}\frac{\mathrm{\Phi }_d}{2},$$
(19)
where

$$f(\mathrm{\Delta }m^2,\theta )=\frac{2\mathrm{cos}2\overline{\theta _m^0}\mathrm{sin}^22\theta }{1+\mathrm{cos}2\overline{\theta _m^0}\mathrm{cos}2\theta }.$$
(20)
Here we have taken $`P^0=(1+\mathrm{cos}2\overline{\theta }_m^0\mathrm{cos}2\theta )/2`$ neglecting the Earth matter effect. The dependence of $`f(\mathrm{\Delta }m^2,\theta )`$ on $`\mathrm{tan}^2\theta `$ for $`\mathrm{\Delta }m^2`$ from the allowed LMA region is shown in Fig. 1. According to Fig. 1 the function $`f(\mathrm{\Delta }m^2,\theta )`$ has maximum $`f^{max}0.4`$ at about the best fit value of $`\mathrm{tan}^2\theta `$.

Fig. 1. The dependence of $`f(\mathrm{\Delta }m^2,\theta )`$ on $`\mathrm{tan}^2\theta `$ for $`\mathrm{\Delta }m^2=510^5`$ eV2, $`710^5`$ eV2 and $`10^4`$eV2.

The relative change of the probability can be rewritten as

$$\frac{\mathrm{\Delta }P}{P^0}=A_{max}\mathrm{sin}\left(\mathrm{\Phi }_2+\frac{\mathrm{\Phi }_d}{2}\right)\mathrm{sin}\frac{\mathrm{\Phi }_d}{2},$$
(21)
where

$$A_{max}\frac{\mathrm{\Delta }P_e}{P_e}|_{max}=ϵ\xi f(\mathrm{\Delta }m^2,\theta )$$
(22)
is the maximal value of $`\mathrm{\Delta }P/P^0`$ which corresponds to both sines being 1. The effect increases with decrease of $`\mathrm{\Delta }m^2`$. For minimal allowed value $`\mathrm{\Delta }m^2=510^5`$eV2, and $`\mathrm{tan}^2\theta =0.3`$, $`\xi =0.64`$ we get

$$A_{max}=0.1\%.$$
(23)
In Fig. 2 we show $`A_{max}`$ as the function of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$. The effect weakly depends on $`\theta `$ in the allowed LMA region.

Fig. 2. Maximum value of the relative effect, $`A_{max}(\mathrm{\Delta }P/P)_{max}`$, as the function of $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$.

If detector is placed underwater, neutrinos will cross an additional layer of water, $`l_W`$,which can be rather large. Straightforward calculations lead to result, similar to eq.(19) in which

$$\mathrm{\Phi }_2\mathrm{\Phi }_2^{}=\mathrm{\Phi }_2+\mathrm{\Phi }_W,\mathrm{\Phi }_w\frac{\mathrm{\Delta }m^2l_W}{2E},$$

where $`\mathrm{\Phi }_w`$ is the phase acquired by neutrinos in water.

Let us comment on possible effect of the third neutrino. We consider the scheme in which $`\nu _3`$ is separated from two others by the mass gap $`\mathrm{\Delta }m_{23}^2=\mathrm{\Delta }m_{atm}^2=(23)10^3`$ eV2. (So, the scheme explains the atmospheric neutrino data.) We assume that this neutrino has small admixture of the electron neutrino described by $`U_{e3}<0.15`$. The influence of matter on this mixing is determined by $`ϵ_{atm}=ϵ(\mathrm{\Delta }m_{atm}^2)710^5`$. Consequently, variations of the flux due to cavity effect should be smaller than $`710^5`$. In fact, the observable effects are further suppressed by smallness of $`U_{e3}`$ and averaging of oscillations associated with the third neutrino. For $`\mathrm{\Delta }m_{atm}^2`$ the oscillation length is smaller than $`1`$ km. Also interference between the modes of oscillations driven by $`\mathrm{\Delta }m_{atm}^2`$ and solar $`\mathrm{\Delta }m^2`$ produces negligible effect. So, we conclude that the cavity effect is mainly due to oscillations driven by solar $`\mathrm{\Delta }m^2`$.

## 3 Effect of cavity. Exposure time

Let us consider the time dependence of the electron neutrino flux in, e.g., deep underwater detector. For fixed position of the detector, the time dependence is due to rotation of the Earth. The rotation leads to certain change of the trajectory of solar neutrinos inside the Earth, and consequently, to modification of the oscillation effect. In general, the change of signal with time is both due to the cavity effect and the effect of regeneration in $`P^0`$ (see (5)). The problem is to disentangle these two effects.

Let $`t_{exp}`$ be the exposure time - the time during which the trajectory of solar neutrinos crosses the cavity. The effect of cavity consists of an additional change of the neutrino flux during exposure time $`t_{exp}`$ on the top of usual regeneration effect.

According to (5) and (17) the probability in the absence of cavity equals

$$P^0(t)=\mathrm{sin}^2\overline{\theta }_m^0+\mathrm{cos}2\overline{\theta }_m^0\left[\mathrm{cos}^2\theta ϵ(t)\mathrm{sin}^22\theta \mathrm{sin}^2\frac{\mathrm{\Phi }(t)}{2}\right].$$
(24)
In principle, this probability can be predicted with high enough accuracy. Indeed, the average radial density profile  is a subject to local variations at the level $`(57)\%`$ . Therefore one expects $`(57)\%`$ variations of $`ϵ(t)`$. This uncertainty can be important for very small values of the density contrast: $`\xi (57)\%`$ or very small phases in Eq. (19). Thus, one should look for the deviation of the time dependence of signal from the signal expected according to (24). Identification of the cavity effect depends on specific conditions of observation. If the cavity is close enough to the detector, the distance $`l_1`$, and therefore $`L`$, can change significantly during the exposure time, so that the averaging over the phase $`\mathrm{\Phi }`$ will occur. As a result

$$\overline{P^0}=\mathrm{sin}^2\overline{\theta }_m^0+\mathrm{cos}2\overline{\theta }_m^0\left(\mathrm{cos}^2\theta \frac{ϵ}{2}\mathrm{sin}^22\theta \right)$$
(25)
practically does not depend on time, and whole time dependence is due to presence of cavity.

Due to smallness of relative change of the probability, $`\mathrm{\Delta }P/P`$, one needs to collect the integral effect during long time of observation. Therefore, the effect of cavity is determined by the relative change of probability averaged over the exposure time:

$$\frac{\overline{\mathrm{\Delta }P}}{P}=\frac{1}{t_{exp}}A_{max}𝑑t\mathrm{sin}\left(\mathrm{\Phi }_2^{}(t)+\frac{\mathrm{\Phi }_d(t)}{2}\right)\mathrm{sin}\frac{\mathrm{\Phi }_d(t)}{2}.$$
(26)
In general, the cavity has complicated form, so that $`d=d(t)`$ can change with time during exposure interval significantly. Averaging over $`d`$ gives

$$\frac{\overline{\mathrm{\Delta }P}}{P}=\frac{1}{t_{exp}}\frac{A_{max}}{2}𝑑t\mathrm{cos}\mathrm{\Phi }_2^{}(t).$$
(27)
The size of the cavity should not be much smaller than the oscillation length. For the $`\mathrm{\Delta }m^2710^5`$eV2 the oscillation length is 30 km. Therefore maximum effect should be seen for the cavity size $`l_\nu /215`$ km. For smaller cavities the effect will be suppressed according to the Eq. (19).

$`\mathrm{\Delta }P`$ is the periodic function of the phase $`\mathrm{\Phi }_2^{}`$. Averaging over $`\mathrm{\Phi }_2^{}`$ leads to disappearance of the cavity effect $`\overline{\mathrm{\Delta }P}=0`$. Suppose averaging is absent. By moving the detector one can vary distances $`l_2`$ and $`l_W`$. Selecting the distance $`l_2`$ and $`l_W`$ in such a way that $`\mathrm{\Phi }_2^{}+\mathrm{\Phi }_d/2=\pi /2+\pi k`$ ($`k=`$ integer) one can maximize the effect

$$\mathrm{\Delta }P=ϵ\xi \mathrm{sin}^22\theta \mathrm{cos}2\overline{\theta }_m^0\mathrm{sin}\frac{\mathrm{\Phi }_d}{2}.$$
(28)
Correspondingly, there are zones of $`l_2`$ with high and low sensitivity. Moving the detector by $`l_\nu /2`$ one can cover whole space.

The width of the beryllium neutrino line, $`\mathrm{\Delta }E_{Be}1.5`$ keV , is very small. Consequently, no averaging over the energy occurs. Averaging would require the length of trajectory: $`l_1>l_\nu E_{Be}/\mathrm{\Delta }E_{Be}210^4`$ km which is larger than the Earth diameter.

The sign of the effect depends on the sign of density contrast, the distance between the cavity and the detector, $`l_2`$, and the size of the cavity itself. In particular, if the cavity is close enough to the detector, so that $`\mathrm{\Phi }_2\pi `$, and therefore $`\mathrm{\Delta }P\xi \mathrm{sin}^2\mathrm{\Phi }_d/2`$, the $`\nu _e`$ flux (signal) increases, when neutrinos cross the cavity for the positive contrast $`\xi >0`$ $`\rho _d<\rho _2`$ (which is satisfied for oil or water), and it decreases for the negative contrast. Similar situation is for $`\mathrm{\Phi }_22\pi k`$, $`k=1,2,\mathrm{}`$. The sign of the effect changes if the detector approaches or removes from the cavity by the distance $`l_\nu /2`$.

Detecting the cavity from the first position of the detector one determines the line along which the cavity is situated. Observing the cavity from two different positions of detector one can reconstruct location of the cavity unless another cavity is detected from the second position. The latter can be checked by observation from the third position of the detector.

In general, one should take into account finite size of the region in which beryllium neutrinos are produced inside the sun. The diameter of the production region ($`0.2R_{}`$) has the angular size at the Earth $`\theta _{Be}=910^4`$ radian. Therefore, the cavity with transverse size $`h`$ shields the production region completely if

$$l_2\frac{h}{\theta _{Be}}1100\mathrm{km}\left(\frac{h}{1\mathrm{k}\mathrm{m}}\right).$$
(29)
If the cross-section of the cavity, $`S_d`$, is smaller than $`(l_2\theta _{Be})^2`$ an additional suppression factor appears in the size of the relative change of the flux:

$$\frac{\mathrm{\Delta }P}{P^0}\frac{S_d}{(l_2\theta _{Be})^2}.$$
(30)
The integral effect is determined by the exposure time, which in turn depends on (i) the size and shape of the cavity, (ii) its location, (iii) position of detector and (iv) season of the year. For simplicity we will consider the cavity of the parallelepiped form situated in the horizontal layer, with $`h`$ being the thickness in the vertical plane, $`l_w`$ the width in the horizontal plane, and as before, $`d`$ the size of the cavity in horizontal plane along the direction from the Sun to the detector.

Let $`t_{exp}^{day}`$ be the exposure time during the day, that is, during one crossing of the cavity by the neutrino trajectory. If $`\nu `$ trajectory crosses the cavity $`N`$ days during a year the total exposure time during the year equals:

$$t_{exp}t_{exp}^{day}N.$$
(31)
Let us estimate the exposure time for different possible positions of the cavity.

Suppose the cavity is in the equatorial region. If the detector is situated in the same horizontal plane as the cavity, the diurnal exposure time is determined by the cavity thickness:

$$t_{exp}^{day}\frac{t_{day}h}{2\pi l_2}2.3\mathrm{min}\left(\frac{h}{1\mathrm{k}\mathrm{m}}\right)\left(\frac{100\mathrm{k}\mathrm{m}}{l_2}\right),$$
(32)
where $`t_{day}24`$ hours. If the detector is out of the plane or cavity is inclined, the exposure time can be larger. However, at the same time the length of $`\nu `$ trajectory in the cavity will be shorter. Both factors can compensate each other in the integral number events.

The cavity will shield the detector during $`N=365`$ days, so that the total exposure time during the year can be as large as $`14`$ hours.

Let us consider the cavity near the North pole. Here the search will be possible in spring or in autumn when the sun appears at horizon. At this time the sun rises by $`710^3`$ radian/day. Therefore the sun will be shielded by the cavity during $`N=h/(710^3l_2)`$ days. For thickness $`h=1`$ km we get $`N=1.5`$ days. The diurnal exposure is determined by the width $`l_w`$, so that total exposure time equals

$$t_{exp}=\frac{l_wt_{day}N}{2\pi l_2}day1\mathrm{hour}\left(\frac{l_{width}}{20\mathrm{k}\mathrm{m}}\right)\left(\frac{100\mathrm{k}\mathrm{m}}{l_2}\right)\left(\frac{h}{1\mathrm{k}\mathrm{m}}\right),$$
(33)
Let us consider now the search at latitudes 670 which correspond to the polar circle. In this case in December the sun is at horizontal (south) direction for 18 days a year. The diurnal exposure time is the same as in previous case. As a result, $`t_{exp}=13.5`$ hours. Similar consideration holds for observation in June at the south latitude 670.

## 4 Detection

Clearly, detection of so small ( $`<0.1\%`$) effect at low energies and during restricted time interval is a very challenging task. It may require new technological developments. Here we will give simple estimations which help to understand problems of the method.

Let us consider a possibility to detects the effect of the cavity by large underwater detector-submarine.

The background conditions should be as good as in the underground detectors of Beryllium neutrinos . To suppress the background generated by the cosmic rays, the detector should be placed deep underwater $``$ at about 3500 m below sea level. This certainly restricts applications of the method and makes it more difficult technologically. Indeed, it is unlikely that, e.g., oil reservoirs are very deep although such possibility is not excluded . On the other hand, there is not too many places where the sea becomes very deep near the cost.

The background and therefore the depth could be reduced for detectors which use light elements (e.g., helium) and efficient anti-coincidence shielding.

Let us estimate the change of number of events due to presence of cavity and required size of the detector. For definiteness we will consider detection of neutrinos via the $`\nu e^{}`$scattering in scintillator. The rate of events due to the Beryllium neutrino flux without oscillations is :

$$r_0=20\frac{\mathrm{events}}{\mathrm{kt}\mathrm{hour}}.$$
(34)
Then in presence of neutrino conversion the expecting rate is

$$r=\left[P+\frac{1}{5}(1P)\right]r_0,$$
(35)
where $`P`$ is the survival probability of the electron neutrinos. Here we have taken into account that $`\nu _ee^{}`$ cross section is about 5 times larger than $`\nu _\mu e^{}`$ cross section for $`{}_{}{}^{7}Be`$ neutrino scattering with energy transfer (0.25 - 0.664) MeV. Taking typical value $`P^00.6`$ we get

$$r=14\frac{\mathrm{events}}{\mathrm{kt}\mathrm{hour}}.$$
(36)
Total number of events detected during the year equals

$$N=rM_dt_{exp}^{year},$$
(37)
where $`M_d`$ is the mass of detector.

The effect of cavity consists of the increase (or decrease) of the number of events by amount

$$\mathrm{\Delta }N_c=N\overline{\frac{\mathrm{\Delta }P}{P}}\frac{1}{1+1/4P^0}=0.7N\overline{\frac{\mathrm{\Delta }P}{P}}$$
(38)
(see Eq. (26)). Last factor in this equation reflects dumping effect due to contribution of neutral current interactions. This factor is absent for detection via the charged current interactions.

Let us find the number of events needed to establish an excess (or deficit) of the number of event due to presence of the cavity at $`k\sigma `$ level ($`k`$ is integer). We take an ideal situation of 100 % efficiency of detection, absence of background, absence of averaging, etc.. In this case the possibility to identify the effect of cavity will be determined by statistical fluctuation of the neutrino signal during the exposure time: $`\mathrm{\Delta }N_{stat}=\sqrt{N}`$. The effect of cavity should be larger than the fluctuation. The effect of cavity at $`k\sigma `$ level will be achieved if $`\mathrm{\Delta }N_ck\mathrm{\Delta }N_{stat}=k\sqrt{N}`$. From this condition we find using Eq. (38):

$$N2k^2\left(\overline{\frac{\mathrm{\Delta }P}{P}}\right)^2.$$
(39)
For $`k=3`$ ($`3\sigma `$ level) and $`\mathrm{\Delta }P/P=0.001`$ this inequality gives

$$N1.810^7\mathrm{events}.$$
(40)
Consequently, the mass of detector should be

$$M_d\frac{N}{rt_{exp}}=\frac{2k^2}{rt_{exp}}\left(\overline{\frac{\mathrm{\Delta }P}{P}}\right)^2.$$
(41)
Taking the exposure time $`t_{exp}10`$ hours we find $`M=130`$ Mt. In general,

$$M_d=130\mathrm{M}\mathrm{t}\frac{k^2}{9}\frac{10\mathrm{h}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{s}}{t_{exp}}.$$
(42)
For the density of scintillator $`0.8`$ g/cc a detector could have dimensions $``$(500 m)3. (Compare with SuperKamiokande which has diameter 39 m and the height 42 m.) Notice that for such a large detector the problem can appear with fiducial volume. In our estimations of the exposure time we have neglected the finite sizes of the detector. In fact, only part of the so large detector can be shielded by the cavity: that is, the trajectories of neutrinos which cross the cavity will cover only part of the detector during all exposure time. Therefore the effect of cavity will be further suppressed.

We have found that the size of detector (even in the most favorable situation) should be an order of magnitude larger than the size of future underground detectors which are now under consideration . Solution of this problem may require next step of technological developments.

## 5 Discussion and conclusions

We have studied effects of the thin layers of matter on oscillations of solar neutrinos, in particular, on the Be-neutrinos. The effect is proportional to $`ϵ`$, the effective mixing parameter $`\mathrm{cos}2\overline{\theta }_m^0`$ in the production point, and the density contrast.

The oscillations occur in the vacuum dominating regime and the matter effect appears as the small correction to the vacuum oscillations pattern. For the beryllium line the oscillation length is about 20 - 40 km and the matter effect parameter $`ϵ310^3`$. This parameter put the absolute upper bound on the effect for a given energy. The only possibility to further enhance the effect is to have several layers with configurations which add up their effects coherently .

The interesting feature of the oscillation result is that for solar neutrinos it does not depend on the oscillation phase acquired by the neutrinos before the layer with density contrast.

We have considered a possibility to study small scale structures inside the Earth via oscillations of the solar Beryllium neutrinos. Such a study can be sensitive to cavities with density contrast $`O(1)`$ and the size (10 - 20) km situated not too far ($`<`$ hundreds km) from the detector. The cavity can produce a change of the flux by about 0.1 % at most. We have found that exposure time during the year can reach about 10 hours. So, the detection of a cavity at $`3\sigma `$ level would require a detector of about 130 Mt size.

There are some ways to reduce the mass of detector:

1) The exposure time can be increased if the detector moves following the neutrino trajectories which cross the cavity. For this, of course, some a priori knowledge of possible position of the cavity is needed.

Another possibility is to perform observations in the same place during several years.

2) If sterile neutrinos participate in oscillations one can use coherent effects on scattering of neutrinos on nuclei by neutral currents in order to increase the cross section .

Refraction effect increases with neutrino energy: $`ϵE`$. For the boron neutrinos we get $`ϵ0.03`$. Moreover, the cross section of the $`\nu e`$scattering is proportional to neutrino energy. Finally, the background conditions are better. So, one may think to use boron neutrinos instead of beryllium neutrinos. However, the flux of boron neutrinos is 3 orders of magnitude smaller, and therefore total number of events turns out to be one order of magnitude smaller than in the case of beryllium neutrinos. There are other problems: the oscillation length is now about 300 km, and therefore with boron neutrinos one can probe larger structures than with $`Be`$ neutrinos Furthermore, due to continuous spectrum the problem appears with averaging over energies of neutrinos.

Further studies are needed to show feasibility of this method to study of small scale structures inside the Earth and for the geological research.


text:
# Particle Correlations at LEP

## 1 Introduction

The analysis of particle correlations in high energy interactions gives important information about the hadron production mechanism, complementary to studies of global event properties and single–particle distributions.

The LEP $`\mathrm{e}^+\mathrm{e}^{}`$ collider provides an ideal environment for such studies. From 1989 until 1995, LEP operated at centre–of–mass energies around 91 GeV which allowed each of the four experiments, ALEPH, DELPHI, L3 and OPAL, to record more than four million hadronic $`\mathrm{Z}^0`$ decays. From 1996, after the collider energy had been increased above the WW threshold, until the end of LEP four years later, each experiment recorded about ten thousand $`\mathrm{W}^+\mathrm{W}^{}`$ events. A hadronic decay of a Z or W boson leads to some dozen particles in the final state, mostly charged pions and photons from the decay of the $`\pi ^0`$ mesons, but also, to a lesser extent, to kaons, protons and lambda hyperons, which allows to study particle correlations in detail.

Bose–Einstein correlations (BEC) between identical bosons are a well established phenomenon in high energy physics experiments and are often considered to be an equivalent of the Hanbury Brown & Twiss (HBT) effect in astronomy describing the interference of photons emitted incoherently. An alternative approach was proposed by Andersson et. al. taking into acount the dynamics of hadron formation in a coherent production process within the framework of the Lund string model related to the symmetrisation of the quantum–mechanical amplitude.

Bose–Einstein correlations lead to an enhanced production of identical boson pairs with a small four–momentum difference $`\mathrm{Q}^2=(\mathrm{p}_1^\mu \mathrm{p}_2^\mu )^2`$. Traditionally it is studied using a two–particle correlation function $`\mathrm{C}(\mathrm{p}_1,\mathrm{p}_2)=\rho _2(\mathrm{p}_1,\mathrm{p}_2)/\stackrel{~}{\rho }_2(\mathrm{p}_1,\mathrm{p}_2)`$, where $`\rho _2`$ and $`\stackrel{~}{\rho }_2`$ are the two–particle densities with and without BEC, respectively. For the construction of the reference sample $`\stackrel{~}{\rho }_2`$ frequently a MC model without BEC is used. Following the pioneering analysis of Goldhaber, Goldhaber, Lee and Pais (GGLP), a correlation function of type $`\mathrm{C}(\mathrm{Q})=1+\lambda \mathrm{exp}(\mathrm{Q}^2\mathrm{r}^2)`$ is often used to yield a value for r, which is interpreted to be the emitter radius. The factor $`\lambda `$ measures the strength of the BEC effect but sometimes absorbs experimental inpurities. Extra terms in the parametrisation are occasionally used to account for imperfections in the description of other correlations in the reference sample.

## 2 Particle Correlations in $`𝐙^\mathrm{𝟎}`$ Decays

### 2.1 Bose–Einstein Correlations in Pion Pairs

Fig. 1 shows recent L3 measurements of the correlation function for charged (left) and neutral (right) pion pairs using a MC without BEC as the reference samples. For both the charged and neutral pion pairs an enhancement at low Q is clearly visible. Using a parametrisation à la GGLP, the obtained source radius for neutral pions is tending to be smaller than r for charged pions, as qualitatively expected in the Lund string model.

### 2.2 Elongation of the Pion Source

A difference in the Bose–Einstein correlation length longitudinally and transversely with respect to the jet axis in $`\mathrm{e}^+\mathrm{e}^{}`$ annihilation, arises naturally in a model for Bose–Einstein correlations based on the Lund string model.

In the Longitudinal Centre–of–Mass System, which is defined for each pair of particles as the system in which the sum of the two particles’ momenta is perpendicular to the thrust axis of the process, the three–momentum difference $`\stackrel{}{\mathrm{Q}}`$ is decomposed in three components as illustrated in Fig. 2. An analysis with a GGPL parametrisation, separately for the components $`\mathrm{Q}_{\mathrm{long}}`$, $`\mathrm{Q}_{\mathrm{t},\mathrm{side}}`$ and $`\mathrm{Q}_{\mathrm{t},\mathrm{out}}`$, gives access to the transverse and longitudinal emitter radius. In spite of different methods for the construction of the reference sample, the results of DELPHI, L3 and OPAL consistently demonstrate that the particle emission source is elongated along the direction of motion with a ratio of transverse to longitudinal radius in the range 0.6—0.8.

### 2.3 Bose–Einstein Correlations in Kaon Pairs

Bose–Einstein correlations have also been established in charged and neutral kaon pairs at LEP. The left–hand side of Fig. 3 shows the correlation function of selected pairs of charged kaons obtained by OPAL. The kaons were identified using information of the specific ionisation energy loss (dE/dx) in the large volume jet chamber. For neutral kaons an enhanced production at low Q is expected if the C parity of the system is +1 as for two $`\mathrm{K}_\mathrm{S}`$ mesons. On the right–hand side of Fig. 3 the correlation function for pairs of $`\mathrm{K}_\mathrm{S}`$ is shown as obtained by ALEPH. The $`\mathrm{K}_\mathrm{S}`$ mesons were identified by using the $`\mathrm{M}_{\pi ^+\pi ^{}}`$ invariant mass spectrum for candidates with a secondary vertex.

### 2.4 Fermi–Dirac Correlations

Recently it has been proposed to extract an emitter dimension for pairs of equal baryons by utilising the Fermi–Dirac exclusion principle. The correlation function can be parametrised by an equation similar to the GGLP parametrisation with the plus sign replaced by a minus sign. Antisymmetrising the total wave function yields four states, three of which are antisymmetric in space and symmetric in spin. Thus, for an incoherent source, C(Q) should decrease to a value 1/2 in the limit as $`\mathrm{Q}0`$. The left–hand side of Fig. 4 shows the ALEPH results for the correlation function from lambda and anti–lambda pairs. The hyperons were identified with a similar method as described for the $`\mathrm{K}_\mathrm{S}`$ mesons. Using three different methods for the construction of the reference sample (A, B, C), a depletion at low Q is observed in all cases. As shown on the right–hand side of Fig. 4, a depletion for baryons at low Q is confirmed by preliminary results from the OPAL collaboration using pairs of anti–protons, identified by dE/dx information.

### 2.5 Mass Dependence of the Emitter Radius

On the left–hand side of Fig. 5, the measured source radii for pairs of pions (averaged by Alexander et. al.), kaons, lambda hyperons and (anti–)protons are plotted as a function of the hadron mass. A mass hierachy $`\mathrm{r}_\pi >\mathrm{r}_\mathrm{K}>\mathrm{r}_{\mathrm{p},\mathrm{\Lambda }}`$ is visible, although the statement is only justified because of the baryon measurements. The Lund string model in its basic form expects $`\mathrm{r}(\mathrm{m})`$ to increase with m, thus the very small value obtained for $`\mathrm{r}_\mathrm{\Lambda }`$ and $`\mathrm{r}_\mathrm{p}`$ poses a challenge to this model. It has been shown that by applying the Heisenberg uncertainty principles one can derive an expression for $`\mathrm{r}(\mathrm{m})`$ that decreases with the hadron mass m, namely $`r(m)=c\sqrt{\mathrm{}\mathrm{\Delta }t}/\sqrt{m}.`$ In Fig. 5 the prediction is shown with $`\mathrm{\Delta }`$t set to (0.5/1.0/1.5)$``$10-24 sec (lower dashed/solid/upper dashed line), to represent a typical scale for strong interactions. The prediction is able to qualitatively describe the data. As pointed out by Alexander, an enormous energy density of the order of 10–100 GeV/fm3 arises for the baryons if the measured r in fact represents the emitter radius. Another interpretation of the data in the framework of the inside–outside cascade model avoids the potential problem with the extrem energy density. A proportionality between the four–momentum of a produced particle and the four–vector describing its space–time position at the freeze–out is commonly accepted in the description of high–energy collisions. Provided that all particles are emitted from a tube of $``$ 1 fm in diameter at a constant proper time of $``$ 1.5 fm, the model can explain the data as shown in the right–hand side of Fig. 5. In this approach the measured $`\mathrm{r}(\mathrm{m})`$ dependence is solely a consequence of the strong correlation between $`\mathrm{x}^\mu `$ and $`\mathrm{p}^\mu `$.

## 3 Inter–WW Bose–Einstein Correlations

In $`\mathrm{W}^+\mathrm{W}^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ events at LEP, the products of the W decays in general have a significant space–time overlap as the separation of the their decay vertices is small compared to characteristic hadronic distance scales. The W boson mass, a fundamental parameter in the Standard Model, is determined from the corresponding jet masses and could potentially be biased if Bose–Einstein correlations between the decay products of the two W bosons exist. A robust framework to test the presence of such inter–WW BEC was proposed in. If the $`\mathrm{W}^{}`$ and $`\mathrm{W}^+`$ decay independently, then $`\mathrm{\Delta }\rho `$(Q)$``$0 for all Q with the test distribution $`\mathrm{\Delta }\rho `$ defined as:

$$\mathrm{\Delta }\rho =\rho ^{WW4q}2\rho ^{W2q}\rho ^{WW_{mix}},$$

with the two–particle densities $`\rho ^{\mathrm{WW}4\mathrm{q}}`$ determined by the $`\mathrm{W}^+\mathrm{W}^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{q}\overline{\mathrm{q}}`$ sample, $`\rho ^{\mathrm{W}2\mathrm{q}}`$ by the hadronic part of semileptonic $`\mathrm{W}^+\mathrm{W}^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{l}\overline{\nu }`$ events and $`\rho ^{\mathrm{WW}_{\mathrm{mix}}}`$ from events build from two independent semileptonic events without the leptonic parts and combining only particles originating from different W’s. In Fig. 6, the L3 collaboration compares the $`\mathrm{\Delta }\rho `$ distribution obtained from the data with two scenarios of the PYTHIA/PYBOEI Monte Carlo model with BEC. In the upper plot the inter–WW BEC in the Monte Carlo model can be seen as an enhancement of like–sign pairs in the low Q region. Unlike–sign pairs (lower plot) are artificially effected by the technical implementation of the inter–WW BEC in PYBOEI. The data are consistent with no inter–WW correlations and the Monte Carlo which describes BEC between particles from different W bosons in the same way as the correlations within the same W is strongly disfavoured. The same is true for similar results from DELPHI as shown on the right–hand side of Fig. 6. Both results are preliminary.

In the Lund model, Bose–Einstein correlations arise when identical bosons are produced close to each other within the same string. Because of the strong correlation of production space–time and momentum of the hadron in the inside–outside cascade, the measured r is interpreted as the distance in the string where the momentum spectra of the particles still overlap. In the absence of colour reconnection effects, particles from different W bosons are not produced in the same string. However, in addition to the coherent correlations inside a string a second correlation effect of an incoherent HBT type could be present. An analysis of the hadron formation within the Lund model shows that the space–time distance of the production vertices for pairs of particles from different strings is of the order of several fm as illustrated in Fig. 7. For such large distances any remaining inter–WW BEC effect would manifest itself only at very low Q values which are hard to exploit with the limited statistics of WW events at LEP.

## 4 Summary

More than 40 years after BEC entered the high energy physics stage, we now have come to a better understanding of the effect. It has become clear that one has to go beyond a naive HBT interpretation by taking into account the dynamics of hadron production. Since no firm theory exists to describe the hadronisation phase, we rely on phenomenological models. To test such models, particle correlations provide us with details complementary to those obtained from global event properties and single–particle distributions.