text:
# Compressibility crossover and quantum opening of a gap for two dimensional disordered clusters with Coulomb repulsion

## 1 Introduction

The interplay of disorder and interactions in two dimensional ($`2d`$) electronic systems is currently a central problem in condensed matter physics loc99 . The importance of interactions is illustrated by several phenomena, notably the discovery of a metallic phase in various two-dimensional devices Kravchenko . Since the metal-insulator transition occurs at a Coulomb energy to Fermi energy ratio $`r_s10`$, Coulomb repulsion should be crucial in understanding this phenomenon. Compressibility measurements indicate that the high $`r_s`$ insulating phase is incompressible comp1  and spatially inhomogeneous comp2 .

In weakly disordered quantum dots, the inverse compressibility gives the spacing between adjacent conductance peaks. The peak spacing statistics, obtained in experiments in which $`r_s1`$, displays fluctuations larger than those predicted from Random Matrix Theory, characterized by a Gaussian-like distribution instead of a Wigner-Dyson distribution Sivan ; Marcus ; Simmel . This suggests a breakdown of the naive single particle picture, assumed in the constant interaction model.

In strongly disordered insulators, the long range nature of interactions leads to a soft Coulomb gap in the single particle density of states at the Fermi level ES , observed in tunneling spectroscopic experiments in $`3d`$ Lee  and $`2d`$ Ashoori  nonmetallic semiconductors. Assuming that single electron hoppings dominate the transport, the Coulomb gap leads to a crossover in the temperature dependence of the resistivity $`\rho (T)`$ from the Mott variable range hopping law ($`\rho (T)=\rho _M\mathrm{exp}(T_0/T)^{1/3}`$ in dimension $`d=2`$) to the Efros-Shklovskii behavior ($`\rho (T)=\rho _{ES}\mathrm{exp}(T_{ES}/T)^{1/2}`$) ES . This crossover has been reported khondaker  not only decreasing temperature but also decreasing carrier density (around $`r_s1.7`$).

In this paper, we numerically investigate the physics of $`2d`$ Coulomb interacting spinless fermions in a random potential for relatively low values of the factor $`r_s`$ ($`3`$), when the single particle spectra display either Wigner-Dyson statistics for weak disorder (diffusion) or Poisson statistics for strong disorder (Anderson localization). As exact diagonalization studies are restricted to small system sizes paper1 , we use a numerical method CIphys1 ; CIphys2  familiar in quantum chemistry as the configuration interaction method (CIM) CIchem . The method consists in diagonalizing the Hamiltonian in an energetically truncated basis built of the low-energy states of the corresponding Hartree-Fock (HF) Hamiltonian. This method allows us to study the ground state and the lowest energy excitations for different system sizes $`L`$ at constant electronic density.

We show that the average inverse compressibility $`\mathrm{\Delta }_2`$ exhibits a smooth crossover from a $`1/L^2`$ towards a $`1/L`$ decay when $`r_s`$ increases from $`0`$ to $`3`$. At the same time, the distribution of $`\mathrm{\Delta }_2`$ evolves from the Wigner-Dyson distribution (or the Poisson distribution, if disorder is strong enough or $`L`$ large enough) towards a Gaussian-like shape. Therefore fluctuations in $`\mathrm{\Delta }_2`$ are determined by Coulomb repulsion rather than from single particle level fluctuations.

We have also studied the energy level spacing between the ground state and the first excitation for strongly disordered clusters: the average indicates a smooth opening of a gap when $`r_s`$ increases, while the distribution exhibits a sharp Poisson-Wigner-like transition at $`r_s=r_s^C1.2`$. The critical threshold $`r_s^C`$ is characterized by a scale invariant gap distribution, reminiscent of the one particle problem Shklovskii  at a mobility edge. However, it is only the distribution of the first spacing which exhibits such a transition, the distributions of the next spacings remain Poissonian and are essentially unchanged when $`r_s`$ varies.

The paper is organized as follows: in Section 2 we introduce the model, in Section 3 we discuss the numerical technique used and its limits of validity; in Section 4 we present our compressibility data; in Section 5 we discuss the statistical properties of the lowest energy excitations; in Section 6 we discuss the implications of our results for hopping conductivity experiments; in Section 7 we present our conclusions.

## 2 The model

We consider a disordered square lattice with $`M=L^2`$ sites occupied by $`N`$ spinless fermions. The Hamiltonian reads

$$H=t\underset{<i,j>}{}c_i^{}c_j+\underset{i}{}v_in_i+U\underset{ij}{}\frac{n_in_j}{2r_{ij}},$$
(1)
where $`c_i^{}`$ ($`c_i`$) creates (destroys) an electron in the site $`i`$, the hopping term $`t`$ between nearest neighbours characterizes the kinetic energy, $`v_i`$ the site potentials taken at random inside the interval $`[W/2,+W/2]`$, $`n_i=c_i^{}c_i`$ is the occupation number at site $`i`$ and $`U`$ measures the strength of the Coulomb repulsion. The boundary conditions are periodic and $`r_{ij}`$ is the inter-particle distance for a $`2d`$ torus. If $`a_B^{}=\mathrm{}^2ϵ/(m^{}e^2)`$, $`m^{}`$, $`ϵ`$, $`a`$ and $`n_s=N/(aL^2)`$ denote respectively the effective Bohr radius, the effective mass, the dielectric constant, the lattice spacing and the carrier density, the factor $`r_s`$ is given by:

$$r_s=\frac{1}{\sqrt{\pi n_s}a_B^{}}=\frac{U}{2t\sqrt{\pi n_e}},$$
(2)
since in our units $`\mathrm{}^2/(2m^{}a^2)t`$, $`e^2/(ϵa)U`$ and $`n_e=N/L^2`$.

## 3 Configuration interaction method

A numerical study of the model (1) via exact diagonalization techniques for sparse matrices is possible only for small systems paper1 , and does not allow us to vary $`L`$ for a constant density. We are obliged to look for an approximate solution of the problem, using the Hartree-Fock orbitals, and to control the validity of the approximations. One starts from the HF Hamiltonian where the two-body part is reduced to an effective single particle Hamiltonian Kato ; Poilblanc ; Schreiber

$$\begin{array}{c}U(\underset{ij}{}\frac{1}{r_{ij}}n_in_j\underset{ij}{}\frac{1}{r_{ij}}c_i^{}c_jc_j^{}c_i),\end{array}$$
(3)
where $`\mathrm{}`$ stands for the expectation value with respect to the HF ground state, which has to be determined self-consistently. For large values of the interaction and large system sizes the single-particle problem (3) is still non-trivial, since the self-consistent iteration can be trapped in metastable states. This limits our study to small $`r_s`$ and forbids us to study by this method charge crystallization discussed in paper1  at a larger $`r_s^W10`$ .

Fig. 1 shows that the HF approximation gives a good estimate of the ground state energy. This approximation becomes worse for (1) larger interactions or (2) smaller disorder. The first characteristic is due to the fact that the ground state HF energy is exact in the limits $`U0`$ and $`U\mathrm{}`$, in which the ground state is given by a single Slater determinant, but deviations from this simple picture are expected at intermediate $`U`$ values. The latter feature can be understood as complicated many-body screening effects, which are effective up to a distance of the order of the localization length, cannot be reproduced within simple HF.

The main advantage of the HF approximation is that it reproduces well the single particle density of statesSchreiber , particularly the soft Coulomb gap at the Fermi energy ES . However, the approximations involved in the HF method are uncontrolled. The mean field HF results can be improved using the configuration interaction method CIphys1 ; CIphys2 ; CIchem . Once a complete orthonormal basis of HF orbitals has been calculated:

$$H_{HF}(|\psi _1,\mathrm{},|\psi _N)|\psi _\alpha =ϵ_\alpha |\psi _\alpha ,$$
(4)
with $`\alpha =1,2,\mathrm{},L^2`$, it is possible to build up a Slater determinants’ basis for the many-body problem which can be truncated to the $`N_H`$ first Slater determinants, ordered by increasing energies. The two-body Hamiltonian can be written as

$$H_{\mathrm{int}}=\frac{1}{2}\underset{\alpha ,\beta ,\gamma ,\delta }{}Q_{\alpha \beta }^{\gamma \delta }d_\alpha ^{}d_\beta ^{}d_\delta d_\gamma ,$$
(5)
with

$$Q_{\alpha \beta }^{\gamma \delta }=U\underset{ij}{}\frac{\psi _\alpha (i)\psi _\beta (j)\psi _\gamma (i)\psi _\delta (j)}{r_{ij}}$$
(6)
and $`d_\alpha ^{}=_j\psi _\alpha (j)c_j^{}|0`$. One gets the residual interaction subtracting Eq. 3 from Eq. 5. This keeps the two-body nature of the Coulomb interaction, and if $`N2`$ it is still possible to take advantage of the sparsity of the matrix and to diagonalize it via the Lanczos algorithm.

Fig. 2 compares HF and CIM results. Labelling the $`N`$-body levels $`E_i`$ ($`i=0,1,2,\mathrm{}`$) by increasing energy we have studied the first spacing $`\mathrm{\Delta }E_0=E_1E_0`$. When $`r_s>1`$, the residual interaction significantly reduces the mean gap, and slightly changes its distribution. Within HF approximation, the first excitation is a particle-hole excitation starting from the ground state: The mean gap reduction is due to correlation effects beyond the particle-hole interaction. For example, the electron (hole) polarizes its neighborhood as it creates fluctuations in the local charge density.

We have checked that CIM results agree with the results given from exact diagonalization with an accuracy of the order $`2\%`$ when one takes into account the $`N_H=10^3`$ lowest energy Slater determinants when $`N=4`$, $`L=8`$, $`W=15`$, and $`r_s=5`$. This means that a basis spanning only $`0.2\%`$ of the total Hilbert space is sufficient for studying the first excitations. For larger $`L`$, exact diagonalization is no longer possible, but one can look the variation of the results when $`N_H`$ increases. In the worst case considered ($`N=16`$, $`L=16`$, $`W=15`$, $`r_s=2.8`$) the accuracy in the first four spacings can be estimated as better than $`10\%`$ when $`N_H=2\times 10^3`$ (see Fig. 3).

Therefore the CIM method allows to study low energy level statistics for $`r_s<3`$. However, its accuracy is not sufficient to determine more sensitive quantities, like a small change of the ground state energy when the boundary conditions are twisted (i.e. the persistent currents).

## 4 Electronic compressibilty

In quantum dots experiments in the Coulomb blockade regime the spacing between consecutive conductance peaks is given by

$$\mathrm{\Delta }_2(N)=E_0(N+1)2E_0(N)+E_0(N1),$$
(7)
with $`E_0(N)`$ $`N`$-body ground state. This quantity is the discretized second derivative of the ground state energy with respect to the number of particles, i.e. the inverse compressibility. In the constant interaction model, which ignores fluctuations in the Coulomb energy

$$Q_{\alpha \beta }^{\gamma \delta }=\frac{e^2}{C}\delta _{\alpha \gamma }\delta _{\beta \delta },$$
(8)
where $`C`$ is the capacitance of the dot. This gives $`E_0(N)=e^2N(N1)/2C+_{k=1}^N\eta _k`$, where $`\eta _k`$ is the $`k`$-th single particle energy at $`U=0`$. Substituting this expression into Eq. (7) one finds $`\mathrm{\Delta }_2=e^2/C+\eta _{N+1}\eta _N`$. Random Matrix Theory describes the single particle level spacings when the disorder is weak. The peak spacing distribution is expected to follow the Wigner-Dyson distribution with fluctuations in $`\mathrm{\Delta }_2`$ whose standard deviation is given by $`\delta \mathrm{\Delta }_2\sqrt{<\mathrm{\Delta }_2^2><\mathrm{\Delta }_2>^2}=\sqrt{4/\pi 1}<\mathrm{\Delta }>0.52<\mathrm{\Delta }>`$, where $`<\mathrm{\Delta }>`$ is the single particle mean level spacing. However, experiments Sivan ; Marcus ; Simmel  performed at $`r_s1`$ found a distribution which is Gaussian-like and has a larger width, up to $`7.5\mathrm{\Delta }`$ Simmel . This suggests that fluctuations in $`\mathrm{\Delta }_2`$ are dominated by electron-electron interactions. Several observed features of the peak spacing distributions have been reproduced using exact numerical diagonalization of the Hamiltonian model (1) Sivan , HF calculations Levit ; Walker ; Berkovits , and a random matrix model for interacting fermionic systems Alhassid .

In this paper we study the inverse compressibility at different system sizes, for a constant filling factor $`n_e=1/9`$. We consider $`N=4,9,16`$ particles on square lattices of size $`L=6,9,12`$ respectively. A Fermi golden rule approximation for the elastic scattering time gives Walker , for $`n_e1`$, $`k_Fl192\pi n_e(t/W)^2`$, $`k_F`$ and $`l`$ denoting the Fermi wave vector and the elastic mean free path, respectively. We consider $`W=5,8,15`$, corresponding to $`k_Fl2.7,1,0.3`$; the first case corresponds to a diffusive system ($`L>l/a2.3`$, with $`a`$ lattice spacing), the last one to a strongly localized system.

The following inverse compressibility data are obtained with the CIM. We have checked that the residual interaction does not change qualitatively HF results. This is consistent with Fig. 1: HF approximation gives a good estimate of the ground state energy and the inverse compressibility is a physical observable which depends only on the ground state energies at different number of particles. Neither a precise knowledge of the ground state wavefuction (as in the calculation of persistent currents) nor excited states energies (as in studies of spectral statistics) are required.

Fig. 4 shows the $`L`$-dependence of the average inverse compressibility, which is well fitted by the power law $`<\mathrm{\Delta }_2(r_s)>L^{\alpha (r_s)}`$, with $`\alpha (r_s)`$ going from $`2`$ to $`1`$ when $`r_s`$ goes from $`0`$ to $`3`$ approximately. The value $`\alpha =2`$ is expected without interaction ($`<\mathrm{\Delta }_2>=<\mathrm{\Delta }>1/L^2`$). The exponent $`\alpha =1`$ can be understood as follows: for $`r_s>1`$, $`\mathrm{\Delta }_2`$ is dominated by the charging energy and therefore $`<\mathrm{\Delta }_2>1/C1/L`$. In the strongly localized regime ($`W=15`$) the $`1/L`$ law is justified due to the Coulomb gap in the single particle density of states ES : According to Koopmans’ theorem note , assuming that the other charges are not reorganized by the addition of an extra charge,

$$\mathrm{\Delta }_2ϵ_{N+1}ϵ_N\frac{1}{L}$$
(9)
due to the Coulomb gap, where $`ϵ_k`$ is the energy of the $`k`$-th HF orbital at a fixed number $`N`$ of particles.

Fig. 5 shows the distribution of inverse compressibilities, for $`L=12`$, $`W=5`$ and $`W=15`$ (inset). At $`r_s=0`$, $`\mathrm{\Delta }_2`$ distributions are close to the Wigner surmise $`P_W(s=\mathrm{\Delta }_2/<\mathrm{\Delta }_2>)=(\pi s/2)\mathrm{exp}(\pi s^2/4)`$ for $`W=5`$ and to the Poisson distribution $`P_P(s)=\mathrm{exp}(s)`$ for $`W=15`$. In the latter case, deviations from the Poisson distribution at small $`s`$ values are due to the finite system size. In both cases at $`r_s=2.5`$ distributions show a Gaussian shape. This can be understood within the Koopmans’ theorem, as the HF energies are given by

$$ϵ_k=\psi _k|H_1|\psi _k+\underset{\alpha =1}{\overset{N}{}}\left(Q_{\alpha k}^{\alpha k}Q_{\alpha k}^{k\alpha }\right),$$
(10)
with $`H_1`$ one-body part of the Hamiltonian (1): Due to the small correlations of eigenfunctions in a random potential, one can reasonably invoke the central limit theorem (see Ref. Berkovits  for a more detailed discussion).

Fig. 6 shows the $`r_s`$-dependence of the inverse compressibility standard deviation $`\delta \mathrm{\Delta }_2`$, normalized to the single particle mean level spacing $`<\mathrm{\Delta }>`$. One can see that $`\delta \mathrm{\Delta }_2/<\mathrm{\Delta }>`$ increases slowly for $`r_s<1`$ and then more rapidly. This results are in contrast with random phase approximation estimates BA ; Blanter , giving fluctuations of the order $`\mathrm{\Delta }`$. On the contrary, we have checked that $`\delta \mathrm{\Delta }_2/<\mathrm{\Delta }_2>`$ has a weak $`r_s`$ dependence for $`r_s>1`$, in agreement with findings of Ref. Sivan . This confirms that peak spacing fluctuations are dominated by capacitance fluctuations instead of single particle fluctuations.

We point out that our data are limited to spinless fermions. Interesting interaction-induced spin effects have been recently reported in experiments Ensslin  and theoretically investigated Brouwer ; Baranger ; Jacquod .

## 5 First energy excitation

We have calculated the first $`N`$-body energy levels $`E_i`$, $`(i=0,1,2,\mathrm{})`$ for different sizes $`L`$, with a large disorder to hopping ratio $`W/t=15`$ imposed to have Anderson localization and Poissonian spectral statistics for the one particle levels at $`r_s=0`$ when $`L8`$. We studied $`N=4,9`$, and $`16`$ particles inside clusters of size $`L=8,12`$, and $`16`$ respectively. This corresponds to a constant low filling factor $`n_e=1/16`$. We studied an ensemble of $`10^4`$ disorder configurations.

The first average spacing $`<\mathrm{\Delta }E_0>`$ is given in Fig. 7. It exhibits a power law decay as $`L`$ increases, with an exponent $`\beta `$ given in the inset. One finds for the first spacing that $`\beta `$ linearly decreases from $`d=2`$ to $`1`$ when $`r_s`$ increases from $`0`$ to $`3`$. The next mean spacings $`<\mathrm{\Delta }E_i>=<E_{i+1}E_i>`$ depend more weakly on $`r_s`$, as shown in Fig. 7.

For $`r_s=0`$, the distribution of the first spacing $`s=\mathrm{\Delta }E_0/<\mathrm{\Delta }E_0>`$ becomes closer and closer to the Poisson distribution $`P_P(s)`$ when $`L`$ increases, as it should be for an Anderson insulator. For a larger $`r_s`$, the distribution seems to become close to the Wigner surmise $`P_W(s)`$ characteristic of level repulsion in Random Matrix Theory, as shown in Fig. 8 for $`r_s=2.8`$ and $`L=16`$. To study how this $`P(s)`$ goes from Poisson to a Wigner-like distribution when $`r_s`$ increases, we have calculated a parameter $`\eta `$ which decreases from $`1`$ to $`0`$ when $`P(s)`$ goes from Poisson to Wigner:

$$\eta =\frac{\text{var}(P(s))\text{var}(P_W(s))}{\text{var}(P_P(s))\text{var}(P_W(s))},$$
(11)
where $`\text{var}(P(s))`$ denotes the variance of $`P(s)`$, $`\text{var}(P_P(s))=1`$ and $`\text{var}(P_W(s))=0.273`$. In Fig. 9, one can see that the three curves $`\eta (r_s)`$ characterizing the first spacing for $`L=8,12,16`$ intersect at a critical value $`r_s^C1.2`$. Our data suggest that for $`r_s<r_s^C`$ the distribution tends to Poisson in the thermodynamic limit, while for $`r_s>r_s^C`$ it tends to a Wigner-like behavior glass . At the threshold $`r_s^C`$, there is a size-independent intermediate distribution shown in the inset of Fig. 8, exhibiting level repulsion at small $`s`$ followed by a $`\mathrm{exp}(as)`$ decay at large $`s`$ with $`a1.52`$. This Poisson-Wigner transition characterizes only the first spacing, the distributions of the next spacings being quite different. The inset of Fig. 9 does not show an intersection for the parameter $`\eta `$ calculated with the second spacing. The second excitation is less localized than the first one when $`r_s=0`$, since the one particle localization length weakly increases with energy. It is only for $`L=16`$ that the distribution of the second spacing becomes close to Poisson without interaction, and a weak level repulsion occurs as $`r_s`$ increases.

The observed transition, and the difference between the first spacing and the following ones is mainly an effect of the HF mean field. For the first spacing, the curves $`\eta `$ calculated with the HF data are qualitatively the same. At the mean field level the low energy excitations are particle-hole excitations starting from the ground state. The energy spacing between the first and the second excited states is given by the difference of two particle-hole excitations and a Poisson distribution follows if the low energy particle-hole excitations are uncorrelated.

Notice that the energy of an electron-hole pair is given by $`ϵ_jϵ_iU/r_{ij}`$ and the classical argument for the existence of a gap in the single particle density of states does not apply ES . For example, in the constant interaction model one pays a gap $`e^2/C`$ in the single particle density of states but not in the many-body excitation spectrum ($`\mathrm{\Delta }E_0=\mathrm{\Delta }E_0(U=0)1/L^2`$). Therefore the observed opening of a gap for the first energy excitation is a remarkable phenomenon beyond the predictions of the classical Coulomb gap model.

## 6 Hopping conductivity

In disordered insulators, a crossover ES  in the temperature dependence of the resistivity $`\rho (T)`$ is induced by Coulomb repulsion from the Mott variable range hopping law ($`\rho (T)=\rho _M\mathrm{exp}(T_0/T)^{1/3}`$ in $`2d`$) to the Efros-Shklovskii behavior ($`\rho (T)=\rho _{\mathrm{ES}}\mathrm{exp}(T_{\mathrm{ES}}/T)^{1/2}`$).

The usual argument is to consider the length $`L(T)`$ which maximizes the probability $`P_{ij}`$ of a single particle hop from a state with energy $`ϵ_i`$ to a state with energy $`ϵ_j`$:

$$P_{ij}\mathrm{exp}\left(\frac{2L}{\xi }\frac{ϵ_jϵ_i}{kT}\right),$$
(12)
with $`\xi `$ localization length.

To measure possible delocalization effects, we have calculated for $`W=15`$ the number of occupied sites per particle $`\xi _s=N/_i\rho _i^2`$, where $`\rho _i=\mathrm{\Psi }_0|n_i|\mathrm{\Psi }_0`$ is the charge density of the ground state at the site $`i`$. Around $`r_s1.2`$ and after ensemble average, the maximum increase of $`\xi _s`$ compared to $`r_s=0`$ is negligibly small ($`2\%`$). Therefore noticeable effects come mainly from the mean and the distribution of the single particle density of states around the Fermi level.

At low temperatures, hopping conductivity involves states near the Fermi level. If one takes for the Coulomb gap its average value $`ϵ_{N+1}ϵ_N\mathrm{\Delta }_2(A+Br_s)/L^{\alpha (r_s)}`$ (see Fig. 4), one obtains for the hopping resistivity a smooth and continuous crossover from Mott to Efros-Shklovskii hopping, given by:

$$\rho (T)\mathrm{exp}\left(\frac{T(r_s)}{T}\right)^{1/(\alpha +1)},$$
(13)
where

$$T(r_s)\frac{A+Br_s}{k\xi ^\alpha }.$$
(14)
This prediction neglects the crossover from the Poisson distribution to a Gaussian distribution in the inverse compressibility (see Fig. 5), which could be included by considering for $`\mathrm{\Delta }_2`$ a more typical value than its average, for instance the value $`\mathrm{\Delta }_2(b)`$ for which $`_0^{\mathrm{\Delta }_2(b)}P(\mathrm{\Delta }_2^{^{}})𝑑\mathrm{\Delta }_2^{^{}}=b`$, with $`b=1/2`$ for example. This would introduce a crossover in $`T(r_s)`$ around $`r_s1`$.

The existence of a critical $`r_s`$ value for the opening of the Coulomb gap can be understood similarly to PD99 . The single particle density of states around the Fermi energy $`E_F`$ is given by ES  $`\rho (E)|EE_F|/U^2`$ and the gap size $`\mathrm{\Delta }_g=|E_gE_F|`$ can be estimated from the condition $`\rho (E_g)\overline{\rho }`$, with $`\overline{\rho }1/W`$ mean density of states for $`Wt`$, obtaining $`\mathrm{\Delta }_gU^2/W`$. According to Fermi golden rule, the inverse lifetime of a Slater determinant built from electrons localized at given sites is $`\mathrm{\Gamma }_tt^2(1/W)(N/L^2)`$, with $`N/(WL^2)`$ density of states directly coupled by the hopping term of the Hamiltonian (1). This argument suggests that quantum fluctuations can melt the Coulomb gap for $`\mathrm{\Gamma }_t\mathrm{\Delta }_g`$, giving $`r_s1`$. Therefore a crossover from Efros-Shklovskii to Mott hopping conductivity is expected not only increasing temperature but also increasing carrier density, as observed in khondaker  at $`r_s1.7`$.

Notice that in the above arguments we have considered only the excitations leading to the Efros-Shklovskii behavior for the hopping conductivity, that correspond to the withdrawal of one electron or its addition to the system. However, a single electron hop may reorganize the location of the other particles, inducing complex many particle excitations (see Fig. 2, showing that simple HF is unable to reproduce the mean first many-body energy spacing when $`r_s>1`$). Therefore one cannot exclude that correlated hops involving more than one particle affect the hopping conductivity ES .

## 7 Conclusions

In summary, we have analyzed the inverse compressibility and the first energy excitation statistics for spinless fermions in weakly and strongly disordered squared lattices when $`r_s`$ increases from $`0`$ to $`3`$. On one hand, we have found a crossover in the electronic inverse compressibility from a $`1/L^2`$ towards a $`1/L`$ decay; at the same time its distribution evolves towards a Gaussian shape. On the contrary, our data for the first many-body energy spacing are consistent with a smooth opening of a gap and a sharp interaction-induced Poisson-Wigner-like transition in its distribution at $`r_s^C1.2`$, with a scale invariant distribution at the critical point showing intermediate statistics. For those small values of $`r_s`$, the ground state energy is essentially given by an effective Hartree-Fock mean field. We underline that the intermediate quantum regime discussed in Ref. paper1  and charge crystallization occur for larger $`r_s`$ factors, where the mean-field approximation breaks down and the configuration interaction method ceases to efficiently converge. This does not allow us to study the regime characterizing the experiments reported in Refs. comp1 ; comp2 .

###### Acknowledgements.

Partial support from the TMR network “Phase Coherent Dynamics of Hybrid Nanostructures” of the European Union is gratefully acknowledged.


format: text

text:
# Possible critical point in phase diagrams of interlayer Josephson-vortex systems in high-𝑇_𝑐 superconductors

## Abstract

A critical value in the product of the anisotropy parameter and the magnetic field is observed in interlayer Josephson-vortex systems by extensive Monte Carlo simulations. Below/above this critical value the thermodynamic phase transition between the normal and the superconducting states upon temperature sweeping is first/second order. The origin is the intrinsic pinning effect of the layered structure of high-$`T_c`$ superconductors.

PACS numbers: 74.60.Ge, 74.25.Bt, 74.25.Dw, 74.20.De

In the superconducting state, an external magnetic field applied in the Cu-O plane of a high-$`T_c`$ superconductor induces the so-called Josephson vortices. The center of a Josephson vortex enters into a block layer, the layer between two superconducting Cu-O layers, in order to save the condensation energy of superconductivity . The thermodynamic phase transition and the lattice structure of interlayer Josephson vortices have been attracting considerable interests since the discovery of high-$`T_c`$ superconductivity. Using a London theory, the structure of Josephson-vortex lattice was derived as the compressed hexagons of triangular lattice pointing along the c axis by Ivlev, Kopnin and Pokrovsky . The interlayer shear modulus was shown to be exponentially small, and the shear deformation of a rhombic lattice might arise through a second-order phase transition. However, at higher temperatures fluctuations are more important, and the London theory is generally inaccurate for discussions about phase transitions. Considerable efforts have been made in order to clarify the thermodynamic phase transition in Josephson-vortex systems both experimentally  and theoretically  thereafter. Nevertheless, the understanding for the problem is still not satisfactory yet. Results obtained by different techniques even seem to be contrary to each other. The difficulty in approaching this problem is two-fold. On the experimental side, a small deviation of the direction of the magnetic field from the Cu-O plane can lead to a strong influence of the c-axis component of the magnetic field on thermodynamic properties of the systems, since all the high-$`T_c`$ superconductors are very anisotropic. On the theoretical side, one has to treat simultaneously anisotropic inter-vortex forces, the commensuration between the vortex alignment and the underlying layered structure, and thermal fluctuations.

In the present Letter, we report new results of extensive Monte Carlo (MC) simulations on the thermodynamic phase transition and the lattice structure of interlayer Josephson vortices. Our results suggest the existence of a critical value of product of the anisotropy parameter and the magnetic field, such that below/above this critical value the thermodynamic phase transition between the normal and the superconducting states is first/second order upon temperature sweeping.

The model Hamiltonian used for the present simulations is the so-called 3D anisotropic, frustrated XY model defined on the simple cubic lattice :

$$H=J\underset{x\mathrm{axis}}{}\mathrm{cos}(\phi _i\phi _j)J\underset{y\mathrm{axis}}{}\mathrm{cos}(\phi _i\phi _j)\frac{J}{\gamma ^2}\underset{c\mathrm{axis}}{}\mathrm{cos}(\phi _i\phi _j\frac{2\pi }{\varphi _0}_i^jA_c𝑑r_c).$$
(1)
Here the y axis is along the external magnetic field, and $`𝐲𝐜𝐱`$. The unit length of the simple cubic lattice is the distance $`d`$ between the neighboring Cu-O layers in a cuprate. Therefore, the discreteness in the c axis comes from the underlying layered structure of cuprates, while that in the Cu-O planes is introduced merely for computer simulations. The coupling constant is given by $`J=\varphi _0^2d/16\pi ^3\lambda _{ab}^2`$. The anisotropy parameter is defined by $`\gamma =\lambda _c/\lambda _{ab}`$, and determines the ratio between the couplings in the Cu-O plane and along the c axis. In the present model, fluctuations in amplitudes of superconducting order parameters and in the magnetic induction are neglected.

Details of simulation technique are summarized as follows: The density of flux lines induced by the external magnetic field is $`f=Bd^2/\varphi _0`$. A Landau gauge is adopted so that $`A_x=0`$ and $`A_c=xB`$. The system size is $`L_x\times L_y\times L_c=384d\times 200d\times 20d`$, which is compatible with the filling factor $`f=1/32`$. There are 240 Josephson flux lines in the ground state. Periodic boundary conditions are applied on phase variables in all directions. A typical simulation process is started from a random configuration of the phase variables at a high temperature, such as $`T=1.5J/k_B`$. 30000 and 90000 MC sweeps are used for equilibration and statistics, respectively, at each temperature. The last configuration at a temperature is used as the initial configuration at a slightly lower temperature, where the temperature difference is $`\mathrm{\Delta }T=0.1J/k_B`$. Around the transition temperature, more than one million MC sweeps are adopted at each temperature, and meanwhile the cooling rate is reduced to $`\mathrm{\Delta }T=0.01J/k_B`$. Vortices are identified by counting phase differences around plaquettes.

In order to compare our simulation results with existing experimental observations, we choose to study first a system of anisotropy parameter $`\gamma =8`$, which is near to that of YBa2Cu3O7-δ. The magnetic field corresponding to $`f=1/32`$ in our simulations is much stronger than those in experiments, and we come back to this point later. The temperature dependence of the helicity modulus (a quantity proportional to the superfluid density) along the magnetic field and the specific heat is depicted in Fig. 1. There is a clearly observable $`\delta `$-function like peak in the specific heat at $`T_m0.96J/k_B`$, where the helicity modulus along the direction of magnetic field increase sharply from zero. Shown in the same figure is the temperature dependence of the intensities of Bragg peaks in diffraction patterns at $`𝐪_{xc}^{(1)}=(\pm \pi /8d,0)`$ and $`𝐪_{xc}^{(2)}=(\pm \pi /16d,\pm \pi /d)`$. Therefore, a thermodynamic first-order phase transition occurs at $`T_m`$, where the gauge symmetry and translation symmetry are broken simultaneously, corresponding to the realization of superconductivity and Josephson-vortex lattice respectively.

The lattice structure of Josephson vortices at low temperatures is shown in Fig. 2. The unit cell is rhombic with short axis along the c direction and of a length of $`2d`$, and the long axis along the x direction and of a length of $`32d`$. Josephson vortices are distributed in every block layer for the present parameters $`\gamma =8`$ and $`f=1/32`$. This structure is the same as that predicted by Ivlev, Kopnin and Pokrovsky .

The lattice structure in Fig. 2 is obviously the ground state for $`\gamma 8`$ when the filling factor is fixed at $`f=1/32`$. Therefore we can use it for investigations of thermodynamic properties for large anisotropy parameters by a heating process. The specific heats thus obtained are shown in Fig. 3 for anisotropy parameters $`\gamma =8`$, 9, and 10. The $`\delta `$-function peaks in the curves for $`\gamma =8`$ and 9, is suppressed for $`\gamma =10`$ . In Fig. 4 we display the temperature and anisotropy parameter dependence of the phase difference between nearest neighboring Cu-O layers $`\mathrm{cos}(\phi _n\phi _{n+1})`$. There is a jump in $`\mathrm{cos}(\phi _n\phi _{n+1})`$ for $`\gamma =8`$ and 9, which is smeared out for $`\gamma =10`$. As the jump in $`\mathrm{cos}(\phi _n\phi _{n+1})`$ is nothing but the jump in the Josephson energy in units of $`J/\gamma ^2`$, there exists a latent heat at the transition temperature for $`\gamma =8`$ and 9, but not for $`\gamma =10`$, consistently with the data for the specific heat. The value of the latent heat itself is too tiny, about $`\gamma ^2`$ times smaller than that in $`\mathrm{cos}(\phi _n\phi _{n+1})`$, to be detected directly. On the other hand, from a standard finite-size scaling theory for a first-order phase transition, the height of the $`\delta `$-function like peak in the specific heat is proportional to the system size . Therefore, by using a large system such as the one in our simulations, the $`\delta `$-function like peak in the specific heat becomes observable as in Figs. 2 and 3 for the first-order phase transitions.

We have performed simulations for anisotropy parameters $`\gamma =7,6,\mathrm{}`$, down to the isotropic case of $`\gamma =1`$  fixing the filling factor at $`f=1/32`$, and observed first-order phase transitions for all these anisotropy parameters. Therefore, the present simulation results indicate that there is a critical anisotropy parameter in between $`\gamma =9`$ and $`\gamma =10`$ for $`f=1/32`$, below/above which the phase transition is first/second order.

Now we look for the reason of the suppression of the first-order phase transition when the anisotropy parameter is increased. Suppose a complete commensuration is achieved between the alignment of the Josephson vortices shown in Fig. 2 and the underlying layered structure of high-$`T_c`$ cuprates. In other words, the Cu-O layers do not influence the lattice structure of Josephson vortices, but merely fix its position in the c direction. In such a case, the Josephson-vortex lattice should be rescaled into equilateral triangular lattice using the anisotropy parameter $`\gamma `$, and we have a relation as seen in Fig. 5:

$`(2d)^2=d^2+(d/2f\gamma )^2,`$

which results in

$$f\gamma =\frac{1}{2\sqrt{3}}.$$
(2)
Now we increase the anisotropy parameter from that determined by the above relation when the filling factor $`f`$ is fixed. Since the repulsive force between Josephson vortices in the c direction is reduced, the Josephson-vortex lattice would be compressed in this direction in order to achieve the energy minimum for the new anisotropy parameter. This reconstruction of Josephson-vortex lattice is forbidden by the underlying layered structure of the high-$`T_c`$ superconductor. Therefore, the above relation provides a criterion for onset of the intrinsic pinning effect of the layered structure on the formation of Josephson-vortex lattice. For anisotropy parameters larger than that evaluated by the above relation, the lattice structure of Josephson vortices is determined by both the inter-vortex repulsions and the pinning force of the underlying layered structure. The thermodynamic phase transition associated with the formation of the Josephson-vortex lattice can be different in the two regions divided by the above relation.

Numerically, the critical anisotropy parameter for the filling factor $`f=1/32`$ is evaluated as $`\gamma =16/\sqrt{3}9.24`$ by the relation (2). This estimate coincides well with our simulation results, since first-order phase transitions are observed for $`\gamma 9`$ but not for $`\gamma 10`$. We have also performed simulations for the filling factor $`f=1/25`$, and found the variation of phase transition from first to second order around $`\gamma =8`$. This observation is consistent with the relation (2), since for $`f=1/25`$ one has the critical anisotropy parameter $`\gamma 7.22`$. For $`f=1/36`$, we have observed a first-order phase transition even for $`\gamma =10`$, consistently with the critical value $`\gamma 10.39`$. Namely, our simulation results indicate clearly that the critical anisotropy parameter increases with decreasing filling factor, or magnetic field. Quantitatively, the simple relation (2) seems to give a reasonable estimate on the critical anisotropy parameter.

The same variation of the phase transition should be observed when the anisotropy parameter is fixed while the filling factor, or the strength of the magnetic field, is tuned. The relation (2) can be rewritten as

$$B=\frac{\varphi _0}{2\sqrt{3}\gamma d^2}.$$
(3)
For YBa2Cu3O7-δ with $`\gamma 8`$ and $`d=12\AA `$, the critical magnetic field is estimated as $`B50T`$. Therefore the phase transition in the Josephson-vortex systems in YBa2Cu3O7-δ is first order for magnetic fields available experimentally, according to our present study. For Bi2Sr2CaCu2O8+y with $`\gamma 150`$ and $`d=15\AA `$, the critical magnetic field is evaluated as $`B1.7T`$, which can be checked experimentally.

The phase transition of interlayer Josephson vortices in a layered superconductor has been addressed theoretically by Blatter, Ivlev and Rhyner  and Balents and Nelson . In these theories, the Josephson-vortex lattice melts in-between the layers with the formation of a smectic-like vortex liquid, via a second-order phase transition (see also ). Therefore, for strong magnetic fields the theories give a reasonable scenario for the phenomena in Josephson-vortex systems observed in the present simulations.

The importance of thermal fluctuations in the phase transition should be stressed. For example, about one vortex-antivortex pair is thermally excited per Josephson flux at each xc section at $`T=0.8J/k_B`$, a temperature lower than the corresponding transition point, for $`\gamma =8`$ and $`f=1/32`$ . Energetically, an additional Josephson vortex only costs energy of order of $`J/\gamma ^2`$, which becomes very small when the anisotropy parameter $`\gamma `$ is large . Most of the thermally excited Josephson vortices and antivortices are confined in same block layers and form overhangs in flux lines, or closed loops. As the result, for large anisotropy parameters Josephson flux lines induced by the magnetic field collide with each other in same block layers even below the transition temperature . Therefore, these thermally excited Josephson vortices and antivortices play important roles in smearing out the first-order phase transition in a Josephson-vortex system in a layered superconductor with a large anisotropy parameter .

In summary, from extensive Monte Carlo simulations we have found a critical value in the product of the anisotropy parameter and magnetic field in interlayer Josephson-vortex systems in high-$`T_c`$ superconductors. Below/above this critical value, the thermodynamic phase transition between the normal state and the superconducting state is first/second order. According to the present results, the phase transition in Josephson-vortex systems in YBa2Cu3O7-δ is first order under magnetic fields up to $`B50T`$, while it varies from first to second order in Bi2Sr2CaCu2O8+y as the magnetic field is increased to across $`B1.7T`$.

The authors thank L. Bulaevskii, J. Clem, A. Sudbø and Y. Nonomura for useful discussions. The present simulations are performed on the Numerical Materials Simulator (SX-4) of National Research Institute for Metals (NRIM), Japan.

Fig. 1: Temperature dependence of the helicity modulus along the magnetic field, the specific heat and the intensities $`I_1`$ and $`I_2`$ for the Bragg peaks at $`𝐪_{xc}^{(1)}=(\pm \pi /8d,0)`$ and $`𝐪_{xc}^{(2)}=(\pm \pi /16d,\pm \pi /d)`$ respectively for $`\gamma =8`$ and $`f=1/32`$.

Fig. 2: Josephson vortex lattice for $`\gamma =8`$ and $`f=1/32`$ obtained by MC simulations of a cooling process from a random state at high temperatures.

Fig. 3: Temperature dependence of the specific heat for $`\gamma =8`$, 9 and 10 and $`f=1/32`$. Data for $`\gamma =8`$ and 9 are shifted by constants.

Fig. 4: Temperature dependence of the phase difference between nearest neighboring Cu-O layers for $`f=1/32`$. The lines are for eye-guide.

Fig. 5: Real-space unit cell of the Josephson-vortex lattice in a layered superconductor of a large anisotropy parameter.


text:
# 1 Introduction

## 1 Introduction

The flux of solar neutrinos was first measured in the Homestake mine (see  and references therein) over thirty years ago. Since then, it was realized that the measured flux was significantly suppressed with respect to theoretical predictions. More recently, a handful of different experiments have also succeeded in measuring the solar neutrino flux . All experiments measure a neutrino flux which is significantly suppressed with respect to the theoretical predictions of the most recent version of the Standard Solar Model (SSM) . This thirty year old problem is what is referred to as the “solar neutrino puzzle.”

There are different types of solutions to the solar neutrino puzzle. At first sight, it appears natural to suspect that the SSM predictions for the solar neutrino flux are slightly off, and/or that the experiments have underestimated their systematic effects, given that detailed models of the Sun and neutrino experiments are highly non-trivial. However, SSM independent analyses of the neutrino data (see  for a particularly nice and simple example), together with independent experimental evidence in favour of the SSM , seem to indicate that the above solution to the puzzle is strongly disfavoured.

The best solution to the solar neutrino puzzle involves extending the Standard Model of particle physics by assuming that the neutrinos have mass and that they mix, i.e., neutrino mass eigenstates are different from neutrino weak eigenstates. This possibility has become particularly natural in light of the recent strong evidence for $`\nu _\mu `$ oscillations from the atmospheric neutrino data at SuperKamiokande .

Nonetheless, in order to firmly establish that the solution to the solar neutrino puzzle involves physics beyond the Standard Model, it is necessary to come up with SSM independent, robust experimental evidence for, e.g., solar neutrino oscillations. Indeed, these “Smoking Gun” signatures of solar $`\nu _e\nu _{\mathrm{other}}`$ oscillations are among the present goals of the SuperKamiokande experiment, via the measurement of the day-night asymmetry of the solar neutrino data  and the recoil electron energy spectrum , and the SNO experiment , via the measurement of the charged to neutral current ratio, the day-night asymmetry of the data, and the recoil electron kinetic energy spectrum.

Other goals of this and the next generation of neutrino experiments are, if solar neutrino oscillations are established, to determine neutrino oscillation modes and measure masses and mixing angles. The current data allow for different $`\nu _e`$ oscillation modes and a handful of disconnected regions in the mass–mixing-angle parameter space (see  for two-flavour analyses and also  for an extension to the “dark side” of the parameter space).

Experiments dedicated to measuring the flux of low-energy solar neutrinos ($`E_\nu =O(1001000)`$ keV) are going to be extremely useful, and perhaps crucial, in order to fully solve the solar neutrino puzzle. It was recently shown that future experiments (Borexino  and, perhaps, KamLAND ) dedicated to measuring the flux of 7Be neutrinos (produced by 7Be+$`e^{}^7`$Li$`+\nu _e`$ inside the Sun) should be able to establish or exclude the “just-so” solution  to the solar neutrino puzzle via the study of the seasonal variations of the neutrino flux , and establish or exclude the LOW MSW solution  via the study of the zenith angle dependence of their data . Furthermore, the measurement of a sizable 7Be neutrino flux would significantly disfavour the SMA MSW  solution to the solar neutrino puzzle, especially in the case of $`\nu _e\nu _s`$ oscillations (where $`\nu _s`$ is a sterile neutrino, i.e., a standard model singlet), and significantly constrain the SSM independent analysis, which require the flux of 7Be neutrinos to be virtually absent . Finally, we have shown  that, in the advent that the background rates at Borexino and/or KamLAND are exceptionally low, it should be possible to measure a nonzero component of $`\nu _{\mu ,\tau }`$ in the solar neutrino flux by analysing the recoil electron kinetic energy spectrum.

Another exciting possibility is that of measuring the “fundamental” $`pp`$-neutrinos, which are produced in the interior of the Sun by proton-proton fusion ($`p+p^2`$H$`+e^++\nu _e`$) in a real time experiment. Future experiments, such as HELLAZ, HERON, LENS, etc (see  for an overview) are being designed to do just that. The flux of $`pp`$-neutrinos is particularly constrained by the photon flux i.e., the Sun’s luminosity, which, of course, is very well measured on the Earth. These lowest energy solar neutrinos ($`E_\nu \begin{array}{c}<\hfill \\ \hfill \end{array}420`$ keV) are not only the most abundant ones, but also have the best known flux. Their energy spectrum is also very well known, since it is dictated by the particularly well studied $`p+p`$ nuclear fusion reaction. Among these proposed experiments, HELLAZ  will be able to determine the incoming neutrino energy in an event-by-event basis and have the unique opportunity of studying the solar neutrino spectrum and the recoil electron kinetic energy spectrum separately. Similar to what was shown for 7Be neutrinos , the authors of  showed that HELLAZ may be able to measure a nonzero component of $`\nu _{\mu ,\tau }`$ in the solar $`pp`$-neutrino flux by analysing the recoil electron kinetic energy spectrum independent of the SSM prediction for the solar neutrino flux.

In this paper we extend the analysis done in , and study the flavour composition of the flux of $`pp`$ and 7Be neutrinos using the recoil electron kinetic energy spectrum. In particular we will address the capability of future low-energy solar-neutrino experiments to see evidence for $`\nu _{\mu ,\tau }`$ coming from the Sun, and, in light of such evidence, exclude more “exotic” oscillation scenarios, such as $`\nu _e\nu _s`$ or $`\nu _e\overline{\nu }_{\mathrm{any}}`$ oscillations.

Our presentation is organised as follows: Sec. 2 describes the flavour dependent recoil kinetic energy distribution of events at Borexino and HELLAZ. Sec. 3 presents the technique for determining the presence of $`\nu _{\mu ,\tau }`$ coming from the Sun, independent of the SSM prediction for the neutrino flux. We present simulations for both Borexino and HELLAZ and show how such a determination can be improved once we take the SSM prediction for the neutrino flux into account. Sec. 4 describes how the same procedure can be used to exclude the presence of antineutrinos or sterile neutrinos in the solar neutrino flux. In Sec. 5, we conclude.

## 2 Recoil Electron Kinetic Energy Spectrum

In this section, we discuss the differences in the recoil electron kinetic energy spectra among different neutrino species. Low-energy solar neutrinos are detected via “$`\nu `$$`+e^{}`$$`\nu `$$`+e^{}`$ elastic scattering in the experiments which will be considered here. By “$`\nu `$” in the previous sentence, one actually means any of $`\nu _e`$, $`\nu _\mu `$, $`\nu _\tau `$, $`\overline{\nu }_e`$, $`\overline{\nu }_\mu `$, or $`\overline{\nu }_\tau `$. Because $`\nu _\mu `$ and $`\nu _\tau `$ are indistinguishable as far as the reaction above is concerned, we will refer to both as $`\nu _\mu `$.

The kinetic energy distribution of the recoil electrons, for a given incoming neutrino energy $`E_\nu `$ is very well known and given by

$$\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}=\frac{2G_F^2m_e}{\pi }\left[A^2+B^2\left(1\frac{T}{E_\nu }\right)^2AB\frac{m_eT}{(E_\nu )^2}\right],$$
(2.1)
where $`m_e`$ is the electron mass, $`T`$ is the kinetic energy of the recoil electron, and $`G_F`$ is the Fermi constant. The parameters $`A`$ and $`B`$ are given in Table 1. The sign difference in the term $`1/2`$ is a consequence of the presence (absence) of $`W`$-boson exchange and the interchange of $`A`$ and $`B`$ between neutrino and anti-neutrino cases is a consequence of the “handedness” of the weak interactions. Eq. (2.1) is a tree-level expression, but higher order corrections are known to be very small , especially for the neutrino energies of interest, and will be neglected throughout.

Borexino (under construction) is an ultra-pure liquid scintillator tank which detects the scintillating light produced by the recoil electron absorbed by the medium. For more details see . It is sensitive to recoil electron kinetic energies greater than 250 keV, and is therefore sensitive to the (almost) monochromatic 7Be neutrinos with $`E_\nu =862`$ keV. The expected resolution for the kinetic energy measurement varies from roughly 12%, for $`T=T_{\mathrm{min}}=0.25`$ MeV, to 7% for $`T=T_{\mathrm{max}}=0.66`$ MeV . They expect 53 events/day in the SSM (BP95) together with 19 background events/day with the anticipated radiopurity of the scintillator of $`10^{16}`$g/g for U/Th, $`10^{18}`$g/g for 40K, and 14C$`/^{12}`$C$`=10^{18}`$ and no radon diffusion. It is remarkable, however, that the München group of Borexino achieved a radiopurity for an organic liquid (Phenyl-ortho-xylylethane) better than $`1.0\times 10^{17}`$g/g ; this is an upper bound on the contamination, limited by the sensitivity of the neutron activation measurement and hence the actual radiopurity may be even better. In this paper, we ignore the background to the 7Be solar neutrino signal at Borexino. This is probably an overoptimistic assumption, but could be realised in future upgrades given the above-mentioned achievement.

HELLAZ (proposed) is a large time projection chamber (TPC) filled with roughly 2000 m3 of cool helium gas ($`6`$ tons at 5 atmos, 77 K), which serves as the target for $`\nu `$-$`e`$ scattering. The recoil electron propagates in the gas medium before being absorbed, leaving a track of ionization electrons. These are then collected, yielding information about the kinetic energy and the flight direction of the recoil electron. HELLAZ is sensitive to recoil kinetic energies greater than $`50`$ keV, and can therefore “see” most of the $`pp`$-neutrino spectrum. Most importantly, since not only the recoil kinetic energy of scattered electrons is measured but also their direction, it is possible to reconstruct the incoming neutrino energy, given that the position of the Sun in the sky is known, via the simple kinematic relation

$$T=m_e\frac{2\mathrm{cos}^2\theta }{(1+m_e/E_\nu )^2\mathrm{cos}^2\theta },$$
(2.2)
where $`\theta `$ is the recoil electron scattering angle with respect to the incoming neutrino direction in the laboratory frame. Incidently, from Eq. (2.2) it is very easy to compute the maximum value of the recoil electron kinetic energy, $`T_{\mathrm{max}}=T(\theta =0)=E_\nu /(1+m_e/(2E_\nu ))`$. HELLAZ expects to measure the recoil electron kinetic energy with a resolution which varies roughly from 2% to 4% and the incoming neutrino energy with a resolution which varies between 5% and 12% . They expect around 7 events/day from $`pp`$ neutrinos in the SSM with negligible background. The major sources of background at HELLAZ are radioactive impurities from 232Th and 238U in the structure of the TPC. However, because of the detector’s total event reconstruction capabilities (including directional information), very efficient background rejection schemes are possible (see  and references therein for further information).

The issue we would like to concentrate on is whether the shapes of the recoil electron kinetic energy distributions for different (anti)neutrino species are statistically different at Borexino and HELLAZ. With this in mind, Fig. 1 depicts the normalised distribution of events at HELLAZ (left) and Borexino (right).***In addition to $`pp`$-neutrinos, HELLAZ is also sensitive to 7Be neutrinos, as well as the $`pep`$-neutrinos and the neutrinos coming from the CNO-cycle. 7Be neutrinos can be clearly separated from $`pp`$-neutrinos, while the number of $`pep`$ and CNO-cycle neutrino generated events is expected to be less than 10% that of $`pp`$-neutrinos. Borexino is sensitive to, in addition to 7Be neutrinos (with $`E_\nu =862`$ keV), a fraction of the $`pep`$ and the CNO-cycle neutrinos, which produce approximately 10% as many events as 7Be neutrinos. We assume throughout, for simplicity, that only $`pp`$ (7Be) neutrinos are detected at HELLAZ (Borexino). In the case of Borexino, the data is binned into ten kinetic energy bins, between 250 keV and 650 keV. In the case of HELLAZ, the data is binned into $`4\times 21`$ bins in $`E_\nu \times T`$. The bins have a width of 50 keV in the $`E_\nu `$ direction and central values of 245, 295, 345 and 395 keV, while in the $`T`$ direction they have a width of 10 keV in the range from 50 to 260 keV. The bin sizes have been chosen such that they are roughly the same as the resolution of both detectors. In order to integrate over the incoming neutrino energy at HELLAZ, the (normalised) BP98 $`pp`$-spectrum presented at  was used.

Many important features of the recoil electron kinetic distributions are worthwhile to point out. First of all, it is quite clear that the spectrum produced by $`\overline{\nu }_e`$-$`e`$ is much steeper than all the other ones.Some of these features were pointed out in . Second, the $`\nu _e`$ and $`\nu _\mu `$ generated spectra have opposite slopes when the neutrino energy is small enough, while their shapes start to look more and more similar as the neutrino energy increases. Finally, the spectra produced by $`\nu _\mu `$-$`e`$ and $`\overline{\nu }_\mu `$-$`e`$ scattering are extremely similar, especially at very low energies.

All of these features can be readily understood from Eq. (2.1). First, it is convenient to write the expression for the normalised recoil electron kinetic energy distributions,

$$\frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}y}=N[\frac{A}{B}+\frac{B}{A}(1y)^2\frac{m_e}{E_\nu }y,],$$
(2.3)
where $`y=T/E_\nu `$, and $`N`$ is a normalisation constant, such that $`\frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}T}`$ from $`T_{\mathrm{min}}`$ to $`T_{\mathrm{max}}`$ equals unity.

For $`\nu _e`$, $`\frac{A}{B}3`$, while $`\frac{B}{A}1/3`$. For $`\overline{\nu }_e`$ the situation is reversed, while for both $`\nu _\mu `$ and $`\overline{\nu }_\mu `$, $`\frac{A}{B}\frac{B}{A}1`$.In this case, the normalisation constant $`N`$ is negative. Note, curiously enough, that the reason for the similarity between the $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ cases is simply due to the accidental fact that $`\mathrm{sin}^2\theta _W`$ is close to $`1/4`$.§§§The fact that there is a sign difference between $`g_L`$ and $`g_R`$ for muon-type (anti)neutrinos is irrelevant, since these coefficients either appear as squares or as the product $`g_Lg_R`$. This similarity is even more pronounced at very low energies, when the $`\frac{m_e}{E_\nu }`$ term dominates over the $`\frac{A}{B}`$ and $`\frac{B}{A}`$ terms.

Keeping in mind that $`0.59\begin{array}{c}<\hfill \\ \hfill \end{array}\frac{m_e}{E_\nu }\begin{array}{c}<\hfill \\ \hfill \end{array}2.3`$ in the neutrino energy range of interest and $`0<y(1+m_e/2E_\nu )^1<1`$, one may write approximate expressions

$`\left({\displaystyle \frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}y}}\right)_{\nu _e}3{\displaystyle \frac{m_e}{E_\nu }}y,`$ (2.4)
$`\left({\displaystyle \frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}y}}\right)_{\overline{\nu }_e}3(1y)^2{\displaystyle \frac{m_e}{E_\nu }}y,`$ (2.5)
$`\left({\displaystyle \frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}y}}\right)_{\overline{\nu }_\mu }\left({\displaystyle \frac{\mathrm{d}\overline{\sigma }}{\mathrm{d}y}}\right)_{\nu _\mu }1+(1y)^2+{\displaystyle \frac{m_e}{E_\nu }}y.`$ (2.6)
In the limit $`\frac{m_e}{E_\nu }1`$, all three distributions are quite different (see, e.g., Fig. 1(B) in ). The $`\nu _e`$ case is roughly flat, the $`\overline{\nu }_e`$ case ranges from 3 at $`y=0`$ to 0 at $`y=1`$ and the $`\nu _\mu ,\overline{\nu }_\mu `$ case ranges from $`3/2`$ at $`y=0`$ to $`3/4`$ at $`y=1`$.

For $`\frac{m_e}{E_\nu }\begin{array}{c}>\hfill \\ \hfill \end{array}1`$, things are slightly more complicated, but still easy to understand. For example, the slope of the distributions for small values of $`y`$ are, up to normalisation factors, $`\frac{m_e}{E_\nu }`$, $`\left(6+\frac{m_e}{E_\nu }\right)`$ and $`+\left(\frac{m_e}{E_\nu }2\right)`$ for $`\nu _e`$, $`\overline{\nu }_e`$ and $`\nu _\mu ,\overline{\nu }_\mu `$ respectively. It is then easy to note that the $`\overline{\nu }_e`$ slope is significantly more negative than the other two, and that, in the case of $`\nu _\mu ,\overline{\nu }_\mu `$ the slope is actually positive if $`E_\nu `$ is small enough. This is indeed what one observes in Fig. 1.

As the incoming neutrino energy increases, the distributions generated by $`\nu _e`$ and $`\nu _\mu ,\overline{\nu }_\mu `$ look more and more similar. One hint of this behaviour is that the slope of the $`\nu _e`$ case increases (decreases in absolute value), while the slope of the $`\nu _\mu ,\overline{\nu }_\mu `$ decreases. One can easily estimate that for $`0.9`$ MeV$`\begin{array}{c}<\hfill \\ \hfill \end{array}E_\nu \begin{array}{c}<\hfill \\ \hfill \end{array}1.0`$ MeV the shapes of the $`\nu _e`$ and $`\nu _\mu `$ induced recoil kinetic energy distributions are most similar. Indeed, for 7Be neutrino energies, one can already note that the difference between the $`\nu _e`$ and the $`\nu _\mu `$ cases is similar to the difference between the $`\nu _\mu `$ and the $`\overline{\nu }_\mu `$ cases.

## 3 Measuring a $`\nu _{\mu ,\tau }`$ Component in the Solar Neutrino Flux

In this section, we address the question whether the shapes of the recoil electron kinetic energy distributions presented in Sec. 2 are statistically different at Borexino or HELLAZ. In the affirmative case, there is hope that one may be sensitive to a “contamination” of other neutrino types in the solar neutrino flux by analysing the shape of the recoil kinetic energy spectrum. We consider this an “appearance experiment” of the “wrong” types of neutrinos from the Sun. In this section we will only consider the case of $`\nu _e\nu _\mu `$ oscillations.

In the advent of neutrino oscillations, a mixture of different neutrino weak eigenstates reaches the Earth. Given an electron-type neutrino survival probability $`P_{ee}`$, a fraction $`P_{ee}`$ of all the neutrinos arriving at the detector are $`\nu _e`$, while a fraction $`1P_{ee}`$ are $`\nu _\mu `$. The recoil electron kinetic energy distribution will, therefore, be given by

$$\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}=P_{ee}\times \left(\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}\right)_{\nu _e}+(1P_{ee})\times \left(\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}\right)_{\nu _\mu }.$$
(3.7)
Note that, in general, $`P_{ee}`$ is a function of the neutrino oscillation parameters (the mass-squared differences of the neutrino mass eigenstates and the neutrino mixing angles) and the neutrino energy.

We simulate “data” at Borexino and HELLAZ for different values of $`P_{ee}`$. We use the distributions presented in Sec. 2, while the flux of $`pp`$ and 7Be neutrinos are taken from the SSM . In the case of Borexino, the energy dependence of $`P_{ee}`$ is irrelevant, given the monochromatic nature of 7Be neutrinos. In the case of HELLAZ, we assume that $`P_{ee}`$ is constant inside each individual neutrino energy bin. Following the central idea presented in , we perform a $`\chi ^2`$ fit to the “data” using a linear combination of $`\nu _e`$-$`e`$ scattering and $`\nu _\mu `$-$`e`$ scattering with arbitrary coefficients,

$$C_e\times \left(\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}\right)_{\nu _e}+C_\mu \times \left(\frac{\mathrm{d}\sigma (T,E_\nu )}{\mathrm{d}T}\right)_{\nu _\mu },$$
(3.8)
i.e., we perform a two parameter ($`C_e`$ and $`C_\mu `$) fit to the data. This measurement procedure is independent of the SSM prediction for the neutrino flux. Therefore, if a nonzero coefficient of the $`\nu _\mu `$-$`e`$ scattering distribution is measured, one can claim to have detected evidence for neutrinos other than $`\nu _e`$ coming from the Sun. This “appearance” result certainly qualifies as a smoking gun signature for neutrino oscillations.

Fig. 2 (long, thin error bars) depicts the measured value of $`1P_{ee}=\frac{C_\mu }{C_e+C_\mu }`$ in each of the neutrino energy bins defined in Sec. 2 as a function of the input value of $`P_{ee}`$, for 5 years of simulated HELLAZ data. As was mentioned in the previous paragraph, the relevant information one should obtain from the plot is if the measured value of $`1P_{ee}C_\mu `$ is statistically different from zero.

Fig. 3 (right, long, thin error bars) depicts the measured value of $`1P_{ee}`$ as a function of the input value of $`P_{ee}`$, for two years of Borexino running. This is just a repetition of Fig. 2(A) in .***In , a different variable, $`P1P_{ee}`$, was used. Both results are, of course, equivalent. Fig. 3 (left, long, thin error bars) depicts the result obtained at HELLAZ if all energy bins are used in the “data” analysis. This result is only meaningful if $`P_{ee}`$ is roughly constant for neutrino energies ranging from from 220 keV to 420 keV. This happens to be the case for most of the currently preferred regions of the two-neutrino oscillation parameter space, especially LMA, LOW and VAC solutions (see, e.g., ).For the SMA solution, there is a sharp drop in $`P_{ee}`$ at $`E_\nu 0.4`$ MeV, and the three lower bins can be combined without any problem. At HELLAZ, this will show up in the data, as the $`E_\nu `$ spectrum differs from the expected $`pp`$-neutrino spectrum, and hence is not a concern. Clearly, the significance of the measurement is better than the one obtained for individual energy bins (Fig. 2).

Next, the same analysis as above is repeated, except that the SSM prediction for the solar neutrino flux is included in the $`\chi ^2`$ analysis. An uncertainty of 20% (5%) was assumed for the 7Be ($`pp`$) neutrino flux. The theoretical error was considered Gaussian for simplicity.This procedure follows the one used in . The readers are referred to this article for details. Note that the uncertainties assumed here are inflated with respect to the ones quoted in the SSM calculations  (9% and 1%, respectively), in order to render the results very conservative. Since the rates are very high both at Borexino and HELLAZ, the error bars are dominated by the uncertainties in the fluxes and hence they can be approximately scaled according to the assigned flux uncertainties.

The results are presented in Figs. 2 and 3 (short, thick error bars). The significance of the measured value of $`1P_{ee}`$ improves significantly, especially at HELLAZ, because of the small assigned uncertainty on the $`pp`$-neutrino flux. After five years of HELLAZ running, for example, one should be be able to determine a 1-sigma-away-from-zero $`\nu _\mu `$ component in the $`pp`$-neutrino flux even for $`P_{ee}0.9`$. It is also noteworthy that in the case of the SMA MSW solution to the solar neutrino puzzle $`P_{ee}0`$ for 7Be neutrinos, in which case a 4-sigma-away-from-zero evidence for $`\nu _\mu `$ in the solar neutrino flux can be established in only two years of Borexino running! It is important to emphasise that using SSM predictions for the solar neutrino flux is a reasonable thing to do, especially for $`pp`$-neutrinos. As mentioned before, the flux of $`pp`$-neutrinos is very well known because it is tightly related to the flux of light coming from the Sun. It is, therefore, the neutrino flux which is least sensitive to detailed modelling of the Sun’s innards.

Some comments are in order. First, only statistical uncertainties were considered, and there are no background events in our “data.” As discussed in Sec. 2, the assumption of a negligible background rate seems less than realistic at Borexino, but may be possible in future upgrades. It may, however, be a fair assumption in the case of HELLAZ. If the real experimental data contains a sizable number of background events, it is necessary to either subtract the background in a bin-by-bin basis or to somehow model the recoil kinetic energy distribution produced by background events. Analysing either of these procedures, however, is beyond the scope of this paper.

Second, the analysis which does not include the SSM flux predictions is completely model-independent (the only assumption being the electron recoil spectrum as predicted by the standard electroweak theory), while the one which includes the SSM flux predictions is model-dependent. Obviously, one obtains a much better determination of $`P_{ee}`$ with the additional input of the SSM flux predictions. For establishing the “wrong neutrino component” in the solar neutrino flux as a smoking gun signature of the solar neutrino oscillations, the former approach is desired. However, for the purpose of determining the oscillation parameters, the energy dependence of the survival probability, and excluding other neutrino oscillation modes, such as $`\nu _e\nu _s`$ or $`\nu _e\overline{\nu }_{e,\mu ,\tau }`$ (as will be discussed in Sec. 4), it is reasonable to include the SSM predictions in the analysis.

Finally, we point out that the results we obtained for HELLAZ are similar to the ones obtained by J. Séguinot et al . Indeed, we chose neutrino energy bins at HELLAZ which coincide with the ones used in . They also perform two different analyses of their simulated data, one which is independent of the SSM prediction for the solar neutrino flux, and one which assumes the SSM prediction for the flux. However, their data analysis procedure is somewhat different, and they do not take the theoretical uncertainty of the solar neutrino flux prediction into account.

## 4 Testing for the $`\nu _e\nu _s`$ or $`\overline{\nu }_{e,\mu ,\tau }`$ Hypotheses

Although it is most natural to assume that electron-type neutrinos oscillate into some linear combination of muon-type and tau-type neutrinos, there is a logical possibility that electron-type neutrinos might oscillate into standard model singlet sterile neutrinos , or, perhaps, into antineutrinos of all flavours***The original neutrino oscillation paper by Bruno Pontecorvo  did, after all, consider $`\nu _e\overline{\nu }_e`$! (see  and references therein). In this section, we will address the issue of excluding these solar neutrino oscillation modes if the data collected at Borexino and HELLAZ are consistent with $`\nu _e\nu _\mu `$ oscillations.

One can already address these “exotic” oscillation modes with the current experimental data. The flux of electron-type anti-neutrinos from the Sun is particularly constrained by the SuperKamiokande and the LSD experiments : the 95% CL SuperKamiokande upper bound on the flux of $`\overline{\nu }_e`$ from the Sun with energies $`\begin{array}{c}>\hfill \\ \hfill \end{array}6.5`$ MeV is $`\mathrm{\Phi }_{\overline{\nu }_e}<1.8\times 10^5`$ cm-2 s-1, or $`\mathrm{\Phi }_{\overline{\nu }_e}/\mathrm{\Phi }_{SSM}^{{}_{}{}^{8}\mathrm{B}}=0.035`$, where $`\mathrm{\Phi }_{SSM}^{{}_{}{}^{8}\mathrm{B}}`$ is the SSM prediction for the 8B neutrino flux. KamLAND, being a dedicated detector for $`\overline{\nu }_e`$, will improve this bound further down to 0.1% of the SSM flux above reactor anti-neutrino energies ($`E_\nu \begin{array}{c}>\hfill \\ \hfill \end{array}8`$ MeV) after one year of running . There are, however, scenarios in which, for energies below the SuperKamiokande threshold, the $`\nu _e\overline{\nu }_e`$ mixing is quite large ( and references therein). Such a possibility can only be addressed by low-energy solar neutrino experiments.

Below, we discuss the exclusion of electron-type neutrino oscillations into sterile neutrinos or into one of the antineutrino types separately at Borexino and HELLAZ experiments.

### 4.1 $`\nu _e\nu _s`$

The $`\nu _e\nu _s`$ oscillation mode is allowed by the analysis of current solar neutrino data , even though in the case of the MSW solutions to the solar neutrino puzzle, only the equivalent of the SMA MSW solution exists at the 99% confidence level . It is curious to note that, in the case of atmospheric neutrinos, the $`\nu _\mu \nu _s`$ hypothesis is currently somewhat disfavoured .

In the case of $`\nu _e\nu _s`$ oscillations, one expects the recoil electron kinetic energy spectrum to be exactly the same as the one generated by $`\nu _e`$-$`e`$ scattering, since $`\nu _s`$ do not interact with electrons. The only effect of the neutrino oscillations would be to suppress the expected number of events, i.e., the hypothesis of $`\nu _e\nu _s`$ oscillations is identical to assuming that the solar neutrino flux is, somehow, suppressed. Therefore, we attempt to fit the “data” simulated according to Eq. (3.7) in Sec. 3 (remember that the “data” is consistent with $`\nu _e\nu _\mu `$ oscillations) to the trial function Eq. (3.8), where the piece which corresponds to $`C_\mu `$ vanishes identically. This is a one parameter $`\chi ^2`$ fit to $`C_e`$. Note that the only discrimination against $`\nu _s`$ is the recoil energy spectrum, because the rate can be always fitted with the free parameter $`C_e`$. The inclusion of the SSM flux prediction does not help in excluding the $`\nu _e\nu _s`$ oscillation because the free parameter $`C_e`$ makes the predicted flux irrelevant.One can “discover” $`\nu _s`$ by observing a nearly vanishing rate for 7Be neutrinos. In the case of $`\nu _e\nu _{\mu ,\tau }`$ oscillations, neutral-current scattering guarantees at least (when $`P_{ee}=0`$) 21% of the SSM rate. For this purpose, one should rely on the SSM flux prediction. Note that 7Be neutrinos are predicted to have almost completely oscillated to $`\nu _s`$ for the (only available) SMA solution: $`P_{ee}=0.009_{0.005}^{+0.244}`$ .

Fig. 4 depicts the value of $`\chi ^2`$ obtained when one attempts to fit the “data” to the sterile neutrino hypothesis, for 5 years of HELLAZ running (left) and 2 years of Borexino running (right). In the case of HELLAZ, we have assumed that $`P_{ee}`$ is constant over the entire $`pp`$-neutrino energy range. The value of $`\chi ^2`$ is determined using the philosophy employed in , and should be compared to $`N_{\mathrm{bins}}1`$ ($`N_{\mathrm{bins}}`$ is the number of “data” bins). After 5 years of HELLAZ running, one should be able to exclude sterile neutrinos coming form the Sun at more than 99.9% confidence level (CL) if all electron-type neutrino have turned into muon-type neutrinos ($`P_{ee}=0`$). After 2 years of Borexino running, the sterile neutrino hypothesis is only ruled out, at best, at the 89% CL.The situation does improve, of course, it more events are collected at Borexino. After 5 years of Borexino running, for example, one can exclude sterile neutrinos for $`P_{ee}\begin{array}{c}<\hfill \\ \hfill \end{array}0.1`$ at more than 95% CL. The explanation for this is the fact that the recoil electron kinetic energy spectra are very different when one compares the $`\nu _e`$-$`e`$ and the $`\nu _\mu `$-$`e`$ scattering cases at very low energies, i.e., $`pp`$-neutrinos, and similar at $`O`$(MeV) energies, i.e., 7Be neutrinos, as discussed in Sec. 2. The exclusion CL decreases with increasing $`P_{ee}`$, and sterile neutrinos are excluded after 5 years of HELLAZ running only at the 77% CL for $`P_{ee}=0.4`$.

### 4.2 $`\nu _e\overline{\nu }`$, Model-independent Fit

In the case of $`\nu _e\overline{\nu }_{e,\mu }`$ oscillations, we perform a two parameter fit to the “data” simulated as in Sec. 3 to a linear combination of the $`\nu _e`$-$`e`$ and $`\overline{\nu }_{e,\mu }`$-$`e`$ scattering recoil kinetic energy distributions. Fig. 5 depicts the value of $`\chi ^2`$ obtained when such a fit is performed, for 5 years of Borexino and HELLAZ running. The value of $`\chi ^2`$ is to be compared to $`N_{\mathrm{bins}}2`$ to determine exclusion confidence levels. As advertised in Sec. 2, the $`\nu _\mu `$-$`e`$ and $`\overline{\nu }_\mu `$-$`e`$ scattering cases produce almost identical recoil kinetic energy spectra, and are almost undistiguishable at HELLAZ. At Borexino, however, the difference between $`\nu _\mu `$-$`e`$ and $`\overline{\nu }_\mu `$-$`e`$ scattering is similar to the difference between the $`\nu _\mu `$-$`e`$ and $`\nu _e`$-$`e`$ cases, as mentioned in Sec. 2 (see Fig. 1), and some discrimination seems possible. Furthermore, upon close inspection, one should note that the shape of the distribution produced due to $`\nu _e`$-$`e`$ scattering is more similar to the $`\nu _\mu `$-$`e`$ case than the $`\overline{\nu }_\mu `$-$`e`$. Therefore, any $`\overline{\nu }_\mu `$ component in the trial function makes the value of $`\chi ^2`$ larger, i.e., the minimum of $`\chi ^2`$ is obtained when the coefficient of the $`\overline{\nu }_\mu `$ component is zero.§§§We only allow nonnegative coefficients of the distribution functions in the fits, for obvious reasons. This is exactly what happens in the case of $`\nu _e\overline{\nu }_e`$ oscillations, at both experiments. Any $`\overline{\nu }_e`$ component in the flux makes the agreement between the theoretical function and the “data” worse, and again the best value of $`\chi ^2`$ is obtained when the coefficient of the $`\overline{\nu }_e`$-$`e`$ scattering distribution is zero.

One can see from Fig. 5 that, after 5 years of HELLAZ data, $`\overline{\nu }_e`$ coming from the Sun can be ruled out at more than 95% CL if $`P_{ee}\begin{array}{c}<\hfill \\ \hfill \end{array}0.2`$, while $`\nu _e\overline{\nu }_\mu `$ oscillations are not constrained at all, even for $`P_{ee}=0`$. After 5 years of Borexino data, both $`\nu _e\overline{\nu }_\mu `$ and $`\nu _e\overline{\nu }_e`$ oscillations are ruled out at more than 95% CL if $`P_{ee}\begin{array}{c}<\hfill \\ \hfill \end{array}0.1`$.

Even if the $`\nu _e\overline{\nu }`$ hypothesis cannot be ruled out at some reasonable CL, one may still be able to place upper limits on the flux of anti-neutrinos coming from the Sun. In the case of $`\nu _e\overline{\nu }_e`$ oscillations, it is straight-forward to place upper bounds on the flux of electron-type antineutrinos at both HELLAZ and Borexino. The 95% CL upper bounds on the $`\overline{\nu }_e`$ flux are depicted in Fig. 6. Of course, for $`P_{ee}\begin{array}{c}<\hfill \\ \hfill \end{array}0.2`$ (0.1) at HELLAZ (Borexino) the upper bound on the flux is meaningless, since the hypothesis of $`\overline{\nu }_e`$ is already ruled out at more than 95% CL. Note that the upper bounds on the antineutrino fluxes are normalised by the SSM prediction for the $`pp`$-neutrino flux $`\mathrm{\Phi }_{SSM}^{pp}=5.94\times 10^{10}`$ cm-2s-1 for the HELLAZ result, and the SSM prediction for the 7Be neutrino flux $`\mathrm{\Phi }_{SSM}^{{}_{}{}^{7}\mathrm{Be}}=4.8\times 10^9`$ cm-2s-1 for the Borexino result. For comparison, the 95% CL SuperKamiokande upper bound is 3.5% of the SSM flux, while KamLAND will improve it to 0.1% after one year of running. However, both of them are only for the 8B neutrinos. Both the HELLAZ and (especially) the Borexino limits obtained from 5 years of “data” are competitive with the SuperKamiokande limit for lower energy neutrinos.

In the case of $`\nu _e\overline{\nu }_\mu `$ oscillations the situation is more ambiguous, especially at HELLAZ.The same is true at Borexino if one assumes the SSM prediction of the total neutrino flux, as will be described later. Not only are the minimum values of $`\chi ^2`$ very small, but in some cases (especially for small values of $`P_{ee}`$) a zero $`\overline{\nu }_\mu `$ flux is ruled out at more than 95% CL. In such cases, it seems that the reasonable thing to do is to measure the antineutrino flux, not determine upper limits! The only exception to this is the case $`P_{ee}=1`$, when the data looks exactly like the SSM prediction, without neutrino oscillations. Indeed, one can not only set upper limits on the antineutrino fluxes, but should also set limits to the $`\nu _\mu `$ flux. Such limits are presented in Table 2.

It is worthwhile to comment that the information contained in Figs. 5, 6, and in Table 2 is also valid for the case of any unknown source of solar antineutrinos of the electron and the muon-types, not only neutrino oscillations. This is because our “data” was analysed assuming that the total flux of solar neutrinos is unknown. We emphasise that $`P_{ee}`$ is the survival probability of electron neutrinos assuming that they oscillate into active neutrinos, i.e., $`\nu _e\nu _\mu `$ oscillations.

### 4.3 $`\nu _e\overline{\nu }`$, SSM-dependent Fit

Next, the same analysis can be repeated assuming that the solar neutrino flux is known within theoretical errors. Again, the value of $`\chi ^2`$ is computedThis procedure follows the one used in . The readers are referred to this article for details. and compared with $`N_{\mathrm{bins}}2+1`$ (the $`2`$ corresponds to the two coefficients that are varied during the minimisation procedure and the $`+1`$ corresponds to the solar neutrino flux constraint). Fig. 7 depicts the minimised values of $`\chi ^2`$ obtained with 5 years of HELLAZ (left) and Borexino (right) “data.” The theoretical uncertainty on the $`pp`$ (7Be) neutrino flux was taken to be 2% (20%); we inflated the theoretical errors by roughly a factor of two from those in BP98.

A comparison between Figs. 5 and 7 reveals that the exclusion confidence levels increase, sometimes significantly. For example, after 5 years of HELLAZ one can exclude $`\nu _e\overline{\nu }_e`$ oscillations for virtually all values of $`P_{ee}`$ at more than 99.9% CL. This is mostly because the $`\overline{\nu }_e`$ has a total cross section which is significantly larger than $`\nu _\mu `$, and the oscillation $`\nu _e\overline{\nu }_e`$ cannot account for the large suppression in the event rate in the “data” (due to $`\nu _e\nu _\mu `$). Even the elusive $`\nu _e\overline{\nu }_\mu `$ case can be excluded at Borexino at more than 95% CL for $`P_{ee}\begin{array}{c}<\hfill \\ \hfill \end{array}0.2`$. Note that at HELLAZ the ability to discriminate between $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ is still quite limited. It is worthwhile to comment that, unlike in the case of model-independent fits in Sec. 4.2, the minimum value of $`\chi ^2`$ is in general obtained for a nonzero coefficient of the $`\overline{\nu }`$-$`e`$ scattering distribution. The reason for this is that, even though the shape of the $`\overline{\nu }`$-$`e`$ scattering recoil electron kinetic energy distribution is “more wrong,” the contribution to the overall cross section is smaller than the $`\nu _e`$-$`e`$ scattering case, and therefore one obtains values of the solar neutrino flux which are closer to the theoretical ones by having a finite $`\overline{\nu }`$ component, decreasing the value of $`\chi ^2`$.

Again, one may set upper limits on the antineutrino flux. As before, there is some ambiguity with regard to setting upper limits for the $`\overline{\nu }_\mu `$ flux, because for almost all values of $`P_{ee}1`$ at both experiments a zero flux is excluded at more than 95% CL. On the other hand, the $`\nu _e\overline{\nu }_e`$ oscillation hypothesis is almost completely ruled out by HELLAZ and the upper limits obtained at Borexino are not much better than the ones depicted in Fig. 6. For this reason, the equivalent of Fig. 6 in the case at hand is not presented.

Table 3 contains the obtained upper limits on the (anti)neutrino fluxes when $`P_{ee}=1`$, i.e., when the data agrees with the predictions of the SSM. Unlike the case of a free total flux analysis, the results presented in Table 3 assume that the total neutrino flux of neutrinos to be detected at HELLAZ and Borexino is the one predicted by the SSM, i.e., there is no “room” for other, yet unknown, low-energy solar neutrino sources. For this reason, of course, the bounds obtained are (in some cases) much more stringent.

Finally, as argued before, we emphasise that fixing the value of the solar neutrino flux to its SSM value is a reasonable thing to do, especially for $`pp`$-neutrinos. In these “exclusion analyses” such a procedure is even more natural, especially if one keeps in mind that a theoretical hypothesis, i.e., $`\nu _e\nu _\mu `$ oscillations plus the SSM computed values for the solar neutrino flux, has been “confirmed experimentally.”

## 5 Conclusions

In order to unambiguously solve the solar neutrino puzzle, and to establish the oscillations of solar neutrinos (if they occur), clear “smoking gun” signatures are required. Such signatures include a large day-night effect, anomalous seasonal variations, or an obvious distortion of the neutrino energy spectrum. Another unambiguous signature is a discrepancy between the number of charged current and neutral current events at SNO, which can be viewed as an “appearance” experiment of $`\nu _{\mu ,\tau }`$. However, SNO can look for this “appearance” signature only for 8B neutrinos with $`E_\nu \begin{array}{c}>\hfill \\ \hfill \end{array}6.5`$ MeV and hence similar studies for lower energy neutrinos such as 7Be and $`pp`$ neutrinos, which are less sensitive to details of the solar model, are important.

We have argued in this paper that a careful analysis of the recoil kinetic energy spectrum at Borexino and HELLAZ serves as another “smoking gun” signature, in the sense that one may be able to infer, independent of the SSM prediction for the solar neutrino flux, the existence of $`\nu _{\mu ,\tau }`$ coming from the Sun. It is worthwhile to emphasise that this is different from distortions in the incoming neutrino energy spectrum. In our case we are describing an “appearance” experiment, while the analysis of the neutrino energy spectrum is a (energy dependent) “disappearance” experiment.

It is important to point out that, in our simple simulations, no background events were included. While this is probably an oversimplification in the case of Borexino, it may well be a good approximation for HELLAZ. Moreover future upgrades of Borexino may reduce the background further according to recent encouraging progress . One should keep in mind that, even if the background rates are significant, the procedure we described may still be useful if the background can be successfully dealt with (one should not underestimate the ability and creativity of experimental physicists!).

We have also included in the analysis the SSM prediction of the flux of solar neutrinos. While the results obtained in this manner are model-dependent (they are not “smoking gun” signatures of neutrino oscillations), we found them very useful. This is a reasonable thing to do especially for $`pp`$-neutrinos, whose flux is constrained well by the solar luminosity. This additional input makes the measurement of the oscillation probability more precise.

Finally, we have argued that, if the data collected at Borexino and HELLAZ is consistent with $`\nu _e\nu _{\mu ,\tau }`$ oscillations, one can try to exclude other neutrino oscillation modes ($`\nu _e\nu _s`$ and $`\nu _e\overline{\nu }_{e,\mu ,\tau }`$) using the same procedure or, at least, to set upper limits on the flux of solar antineutrinos. Again we considered the possibility of constraining the solar neutrino flux to the SSM predicted value. The main result we obtained is that $`\nu _e\overline{\nu }_e`$ oscillations can, in general, be excluded, while the $`\nu _e\overline{\nu }_{\mu ,\tau }`$ case is much more elusive. Nonetheless, Borexino should be able to exclude $`\nu _e\overline{\nu }_{\mu ,\tau }`$ oscillations if the SMA MSW solution to the solar neutrino puzzle happens to be the correct one.

## Acknowledgements

AdG thanks Hiroshi Nunokawa for questioning whether or not the issue of antineutrinos in the solar flux could be addressed by this procedure during a seminar at the Instituto de Física Gleb Wataghin, UNICAMP, Brazil. The work of HM was supported in part by the Director, Office of Science, Office of High Energy and Nuclear Physics, Division of High Energy Physics of the U.S. Department of Energy under Contract DE-AC03-76SF00098 and in part by the National Science Foundation under grant PHY-95-14797 and also by the Alfred P. Sloan Foundation.


text:
# Axial vector current in an electromagnetic field and low-energy neutrino-photon interactions

## I Introduction

The study of neutrino-photon interactions began many years ago when Pontecorvo and Chiu and Morrison  suggested that the process $`\gamma \gamma \nu \overline{\nu }`$ may be important in the analysis of stellar cooling. These interactions are of interest in astrophysics and cosmology.

It is well known that, in the standard model, neutrino-photon interactions appear at the one-loop level. In addition to the heavy gauge bosons which can immediately be integrated out at the desired scales leading to the effective four-fermion interaction, charged fermions (electrons) run in the loop; in particular, the coupling of the photons to these fermions is responsible for the $`\gamma \nu `$ interactions.

A variety of processes of physical relevance in different situations belong to this class of $`\gamma \nu `$ interactions, for which an effective-action description involving only neutrino currents and field strength tensors has partly been found. E.g., Dicus and Repko  have recently derived an effective action for interactions between two neutrinos and three soft photons, which immediately allows for the calculation of all scattering amplitudes of this type. This effective action can be derived via the expectation value of the electromagnetic vector current $`j_\mu =e\overline{\mathrm{\Psi }}(x)\gamma _\mu \mathrm{\Psi }(x)`$; the latter can in turn be calculated very easily from the Heisenberg-Euler effective Lagrangian $`_{\text{HE}}`$ with the aid of the formula $`\overline{\mathrm{\Psi }}(x)\gamma _\mu \mathrm{\Psi }(x)=\delta _{\text{HE}}/e\delta A^\mu (x)`$. In , the effective action of Dicus and Repko was reproduced by a more direct diagrammatic approach.

Incidentally, processes with two photons, for example, neutrino-photon scattering $`\gamma \nu \gamma \nu `$, turn out to be highly suppressed in the standard model  because, according to Yang’s theorem , two photons cannot couple to the $`J=1`$ state. As a result, typical cross sections are exceedingly small and suppressed by factors of $`1/m_W^2`$ (see, e.g., ).

The presence of a medium or external electromagnetic field drastically changes the situation. It induces an effective coupling between photons and neutrinos, even enhancing processes such as $`\nu \nu \gamma `$. Based on the same effective action of Dicus and Repko , it was shown in  that in the presence of an external magnetic field, cross sections for neutrino-photon processes such as $`\gamma \gamma \nu \overline{\nu }`$ and $`\nu \gamma \nu \gamma `$ are amplified by the factor $`\left(m_W/m_e\right)^4\left(B/B_\text{c}\right)^2`$ for soft photon frequencies $`\omega m_e`$ and a weak magnetic field $`BB_\text{c}m^2/e`$ (for extensions to this result, see  and ). The subject of the electromagnetic properties of neutrinos in media was comprehensively studied in  (see also  for neutrinos in magnetized media).

In general, an effective action describing low-energy neutrino-photon processes with an arbitrary number of photons and in the presence of strong external electromagnetic fields appears highly desirable. In other words, we are looking for the analogue of the Heisenberg-Euler effective action which turned out to be extremely useful in QED (see e.g., ).

For this, we start from the effective four-fermion interaction Lagrangian valid at energies very much smaller than the W- and Z-boson masses:

$$_{4f}=\frac{G_F}{\sqrt{2}}\overline{\nu }\gamma ^\mu \left(1+\gamma _5\right)\nu \overline{E}\gamma _\mu \left(g_V+g_A\gamma _5\right)E.$$
(1)
Here, $`E`$ denotes the electron field, $`\gamma _5=i\gamma ^0\gamma ^1\gamma ^2\gamma ^3`$, $`g_V=1\frac{1}{2}\left(14\mathrm{sin}^2\theta _W\right)`$ and $`g_A=1\frac{1}{2}`$ for $`\nu _e`$, where the first terms in $`g_V`$ and $`g_A`$ are the contributions from the W exchange diagram and the second one from the Z exchange diagram. Also, $`g_V=2\mathrm{sin}^2\theta _W\frac{1}{2}`$ and $`g_A=\frac{1}{2}`$ for $`\nu _{\mu ,\tau }`$.

Now concentrating on soft electromagnetic degrees of freedom with energies much smaller than the electron mass, we may integrate out the actual “heaviest” particle, i.e., the electron, and find

$$_{\text{eff}}=\frac{G_F}{\sqrt{2}}\frac{1}{e}\overline{\nu }\gamma ^\mu \left(1+\gamma _5\right)\nu \left(g_\text{V}j^\mu ^A+g_\text{A}j_5^\mu ^A\right),$$
(2)
where the expectation values of the currents are given in terms of the Green’s function in this field: e.g., $`j_5^\mu ^A=\text{i}e\text{tr}\left[\gamma ^\mu \gamma _5G(x,x|A)\right]`$.

Obviously, in order to obtain this effective Lagrangian, one must calculate the expectation values of vector and axial vector currents in a slowly varying electromagnetic field background. As mentioned above, the vector current expectation value can be easily obtained using the well-known Heisenberg-Euler Lagrangian $`_{\text{HE}}`$: $`\overline{E}(x)\gamma _\mu E(x)=\delta _{\text{HE}}/e\delta A^\mu (x)`$. In this way, one obtains a derivative expansion of the vector current around a strong field. For example, the term which is third order in the field and first order in derivatives was used by Dicus and Repko  in their study of $`\nu \gamma \nu \gamma \gamma `$ and cross processes.

With regard to Eq. (2), the derivative expansion of the axial vector current around an arbitrarily strong field is finally required. To our surprise, we could not find such an expression in the vast literature on derivative expansions. Therefore, its derivation remains the final task of our present work.

It is convenient to choose a special gauge condition for the potential $`A_\mu (x)`$ without loss of generality: the so-called Fock–Schwinger gauge***Especially for gradient expansions of the heat kernel or the Heisenberg-Euler Lagrangian, this gauge has proved extremely useful; see, e.g., .

$$(x^\mu y^\mu )A_\mu (x)=0,$$
(3)
which is satisfied by the series

$$A_\mu (x)=\frac{1}{2}(x^\lambda y^\lambda )F_{\lambda \mu }(y)+\frac{1}{3}(x^\lambda y^\lambda )(x^\sigma y^\sigma )_\sigma F_{\lambda \mu }(y)+\mathrm{}$$
(4)
Abbreviating the first term on the right-hand side with $`A_{\text{c}\mu }`$ and the second term with $`a_\mu `$, we may perform a perturbation expansion for the Green’s function with respect to the derivative term $`a_\mu `$:

$$G(x,x^{}|A)=G_c(x,x^{}|A_\text{c})+d^4wG_c(x,w|A_\text{c})ea_\mu (w)\gamma ^\mu G_c(w,x^{}|A_\text{c})+\mathrm{},$$
(5)
where $`G_c(x,x^{}|A_\text{c})`$ is the Green’s function of the electron in the constant field produced by $`A_\text{c}`$. Inserting the expansion (5) into the definition of $`j_5^\mu ^A`$, we obtain, to first order in derivatives,

$$j_5^\mu ^A=\frac{1}{3}_\sigma F_{\alpha \beta }\frac{^2}{k^\sigma k^\alpha }\left[\mathrm{\Pi }_5^{\beta \mu }(k)\right]|_{k=0},$$
(6)
where we encounter the vector–axial-vector amplitude (VA) $`\mathrm{\Pi }_5^{\beta \mu }`$, i.e., the axial analogue of the polarization tensor. This VA amplitude has been calculated very recently by one of the authors in ; an independent confirmation can be found in .

Inserting the presentation of  into Eq. (6), we finally arrive at the first-order gradient expansion of the axial vector current for arbitrarily strong electromagnetic fields valid for slowly varying fields:

$`j_5^\mu ^A`$ $`=`$ $`\text{i}e\text{tr}\left[\gamma ^\mu \gamma _5G(x,x|A)\right]`$ (7)
$`=`$ $`\text{i}{\displaystyle \frac{e^2}{24\pi ^2}}_\sigma F_{\lambda \nu }{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{ds}{s}}e^{ism^2}{\displaystyle \frac{(eas)(ebs)}{\mathrm{sin}(ebs)\mathrm{sinh}(eas)}}`$ (10)
$`\{[C^{\nu \lambda }\left(C^2\right)^{\mu \sigma }+C^{\nu \sigma }\left(C^2\right)^{\mu \lambda }]N_1+[B^{\nu \lambda }\left(B^2\right)^{\mu \sigma }+B^{\nu \sigma }\left(B^2\right)^{\mu \lambda }]N_2`$
$`[3F^{\nu \mu }g^{\lambda \sigma }+F^{\nu \lambda }g^{\mu \sigma }+F^{\nu \sigma }g^{\mu \lambda }]N_3\},`$
where the $`N_i`$ are simple functions of the field strength:

$`N_1={\displaystyle \frac{2\mathrm{sin}(ebs)}{\mathrm{sinh}^2(eas)}}\left[\mathrm{cosh}(eas){\displaystyle \frac{\mathrm{sinh}(eas)}{eas}}\right]+2ebsN_3`$ (11)
$`N_2={\displaystyle \frac{2\mathrm{sinh}(eas)}{\mathrm{sin}^2(ebs)}}\left[\mathrm{cos}(ebs)+{\displaystyle \frac{\mathrm{sin}(ebs)}{ebs}}\right]+2easN_3`$ (12)
$`N_3={\displaystyle \frac{1}{easebs(a^2+b^2)}}\left[a^2{\displaystyle \frac{\mathrm{sin}(ebs)}{\mathrm{sinh}(eas)}}b^2{\displaystyle \frac{\mathrm{sinh}(eas)}{\mathrm{sin}(ebs)}}\right].`$ (13)
For the analysis, we employed a decomposition of the field strength tensor in linearly independent subspaces of electric and magnetic parts:

$$C_{\mu \nu }:=\frac{1}{a^2+b^2}\left(aF_{\mu \nu }+bF_{\mu \nu }^{}\right);B_{\mu \nu }:=\frac{1}{a^2+b^2}\left(bF_{\mu \nu }aF_{\mu \nu }^{}\right),$$
(14)
where the invariants $`a`$ and $`b`$ are defined by

$$a,b=\sqrt{(^2+𝒢^2)^{\frac{1}{2}}\pm };=\frac{1}{4}F_{\mu \nu }F^{\mu \nu },𝒢=\frac{1}{4}F_{\mu \nu }^{}F^{\mu \nu }.$$
(15)
As a consequence, we findWe always use the metric $`g=\mathrm{diag}(+)`$.

$$\left(C^2\right)_{\mu \nu }=\frac{1}{a^2+b^2}\left(F_{\mu \nu }^2+b^2g_{\mu \nu }\right);\left(B^2\right)_{\mu \nu }=\frac{1}{a^2+b^2}\left(F_{\mu \nu }^2a^2g_{\mu \nu }\right).$$
(16)
Concerning our central result for the axial vector current in Eq. (10), we observe its linearity in the derivative of $`F_{\mu \nu }`$ (by construction) which is multiplied by a finite proper-time integral containing the complete dependence on the strength of the fields. We would like to point out that our method of performing the derivative expansion in momentum space (cf. Eq. (6)) is much in the spirit of Ref.  with the essential difference that we expand around an arbitrary value of field strength.

Our representation of $`j_5^\mu `$ in terms of the $`C_{\mu \nu }`$ and $`B_{\mu \nu }`$ tensors is very convenient for performing a weak-field expansion; to second order in $`1/m^2`$, it finally reads:

$`j_5^\mu ^A`$ $`=`$ $`{\displaystyle \frac{e^3}{24\pi ^2m^2}}\left(^\mu 𝒢+(^\alpha F_{\alpha \beta })F^{\beta \mu }\right)`$ (19)
$`+{\displaystyle \frac{e^5}{90\pi ^2m^6}}_\sigma F_{\alpha \beta }[𝒢(F^{\beta \alpha }g^{\mu \sigma }+F^{\beta \sigma }g^{\mu \alpha })+(F^{\beta \alpha }(F^2)^{\mu \sigma }+F^{\beta \sigma }(F^2)^{\mu \alpha })`$
$`(3F^{\mu \beta }g^{\alpha \sigma }+F^{\beta \alpha }g^{\mu \sigma }+F^{\beta \sigma }g^{\mu \alpha })].`$
As a cross-check, one can indeed prove  that the first term can immediately be obtained from the order $`1/m^2`$ term of the famous triangle graph. As a side remark, it should be mentioned that the axial vector anomaly is not present in our result, since we computed the $`1/m^2`$-expansion of $`j_5^\mu `$, whereas the anomaly is mass independent.

Returning to the quest for the general effective action, we can easily insert Eq. (10) into Eq. (2) and a similar expression for the vector current $`j_\mu `$. The well-known latter expression can be derived via the Heisenberg-Euler Lagrangian as mentioned above or, alternatively, via the same method as here proposed for the axial vector current; for this, we would like to stress that Eq. (6) also holds for the vector current with the VA amplitude $`\mathrm{\Pi }_5^{\mu \nu }`$ replaced by the polarization tensor $`\mathrm{\Pi }^{\mu \nu }`$. This tensor in the presence of an arbitrarily constant electromagnetic field was first calculated in  (an appropriate representation can be found in ). Details will be presented elsewhere. We observe that the axial vector current is of even order in the field, while the vector current is of odd order; this is a direct consequence of the Dirac algebra (Furry’s theorem).

To conclude, we have completed the search for a low-energy effective action for neutrino interactions with an arbitrary number of soft photons with wavelengths larger than the Compton wavelength. From a different perspective, we have found the generalization of the Heisenberg-Euler action to the case involving an axial vector coupling.

## Acknowledgments

We would like to thank Professor W. Dittrich for helpful discussions and for carefully reading the manuscript. This work was supported by Deutsche Forschungsgemeinschaft under grant DFG Di 200/5-1.