text:
# One-dimensional non-interacting fermions in harmonic confinement: equilibrium and dynamical properties

## 1 Introduction

Atomic Fermi vapours have been cooled down to degeneracy in several experiments on ultra-cold atoms of 6Li and 40K , using techniques similar to those employed to obtain Bose-Einstein condensation . At very low temperatures collisions between atoms are essentially in the $`s`$-wave channel; however, all even-channel collisions are forbidden for spin-polarized fermions by the Pauli exclusion principle. Since dipolar interactions are also negligible, the Fermi gas can be considered as non-interacting. In the experiments the atoms are trapped by magnetic or optical potentials, providing a confinement which can be well approximated as harmonic. The geometry of the trapping potential can be very anisotropic, thus allowing to reach quasi-1D or quasi-2D situations .

We focus here on the quasi-1D geometry, which is the case of the Paris experiment  and is relevant for experiments on atom lasers and atomic waveguides. The 1D gas of non-interacting fermions serves also as a model for describing the 1D Bose gas with hard-core point-like interactions (Tonks gas, ). Indeed, it was demonstrated by Girardeau that there exists an exact mapping between the bosonic and the fermionc many-body wave functions, $`\mathrm{\Psi }_B=|\mathrm{\Psi }_F|`$, valid also in the time domain . The Tonks gas limit corresponds to the case of total reflection in a two-atom collision , and the interactions in 1D mimic the role of the Fermi statistics.

In this paper we present some results for the exact equilibrium profiles and for the dynamical properties of a 1D Fermi gas in harmonic confinement. In Sec. 2 we illustrate the role of quantum statistics on the particle and kinetic-energy density profiles, and we investigate how the Pauli principle influences the pair correlation function. In Sec. 3 we evaluate the spectrum of the gas both in the collisionless and in the collisional regime. In addition, we relate the result for the frequency of the quadrupole mode to a scaling property of the Hamiltonian used in the long-wavelength limit. Finally, Sec. 4 offers a summary and some concluding remarks.

## 2 Equilibrium properties

We consider a system of $`N`$ non-interacting Fermi particles confined in an external potential which is very tight in two spatial directions. We assume that the gas is in a quasi-1D situation where the transverse degrees of freedom are frozen, that is, the chemical potential is much smaller than the level spacing in the transverse direction. In this case the transverse wave function is the ground state one, which will factorize out from all the following calculations. We can thus focus only on the axial coordinate $`x`$ and define the external confinement along $`x`$ given by a real potential $`V_{ext}(x)`$. For the sake of simplicity we will assume that the eigenfunctions $`\psi _i(x)`$ are real .

The equilibrium properties of the gas can be expressed in terms of the longitudinal one-body density matrix, which at zero temperature reads

$$\rho ^{(l)}(x,x^{})=\underset{i=1}{\overset{N}{}}\psi _i(x)\psi _i(x^{}).$$
(1)
In the following we will give several examples of calculations of the exact spatial profiles obtained as the moments of the density matrix, and we will evaluate also the longitudinal pair distribution function.

### 2.1 Density profiles

The zeroth moment of the one-body density matrix yields the equilibrium density profile $`n(x)`$

$$n(x)=\rho ^{(l)}(x,x^{})|_{x^{}=x}=\underset{i}{}\psi _i|\delta (xx_i)|\psi _i,$$
(2)
where $`x_i`$ are the positions of the fermions and $`|\psi _i`$ are the eigenstates of the Hamiltonian in the coordinate representation. We have shown  that $`n(x)`$ can be written in terms of the Green’s function in coordinate space, $`\widehat{G}(x)=(x\widehat{x}+i\epsilon )^1`$, which is a function of the position operator $`\widehat{x}`$:

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\underset{i=1}{\overset{N}{}}\psi _i|\widehat{G}(x)|\psi _i.$$
(3)
This equation expresses the density profile as the trace of the operator $`\widehat{G}(x)`$ on the first $`N`$ quantum states ($`\mathrm{Tr}_N`$). Using the relation between the partial trace $`\mathrm{Tr}_N`$ of a generic matrix $`Q`$ and the determinant of the inverse matrix $`Q^1`$,

$$\mathrm{Tr}_NQ=[\mathrm{ln}det(Q^1+\lambda 𝕀_N)]/\lambda |_{\lambda =0},$$
(4)
we obtain

$$n(x)=\frac{1}{\pi }\underset{\epsilon 0^+}{lim}\mathrm{Im}\frac{}{\lambda }[\mathrm{ln}det(x+i\epsilon \widehat{x}\lambda 𝕀_N)]_{\lambda =0}.$$
(5)
This alternative expression for the density profile is very practical for performing numerical calculations in the case of harmonic confinement, $`V_{ext}(x)=m\omega _{ho}^2x^2/2`$. For this specific potential the representation for both the position and the momentum operators on the basis $`\psi _i`$ of the eigenstates of the harmonic oscillator is a tridiagonal matrix. Thus the operator $`x+i\epsilon \lambda 𝕀_N`$ is tridiagonal: the determinant of such a matrix can be evaluated recursively through a renormalization procedure .

The numerical results for the particle density profiles are shown in Fig. 1 at various number of fermions, as compared with the predictions of the local-density approximation: the latter neglects the shell structure and the spill-out from the classical border. Because of the symmetry between position $`\widehat{x}`$ and momentum $`\widehat{p}`$ in the harmonic-oscillator Hamiltonian, the momentum distribution has exactly the same behaviour as the particle density illustrated in Fig. 1.

On account of the mapping with the Tonks gas the particle density profile of $`N`$ noninteracting fermions in the 1D harmonic trap is identical to that of $`N`$ bosons with hard-core interactions in the same external potential, while the momentum distribution for such a system has a sharp peak at small wavevectors .

### 2.2 Kinetic energy and momentum flux densities

The first moment and all the other odd moments of the one body density matrix are zero because of the reality of the wave functions. For higher even moments several definitions are avaible in the literature (for a review see e.g. Ref. ); these can be expressed again in terms of $`\widehat{G}(x)`$ . In particular, we shall concentrate on the quantity

$`P_n(x)`$ $`=`$ $`(i)^n\left[{\displaystyle \frac{^n}{x_{}^{}{}_{}{}^{n}}}\rho ^{(l)}(x,x^{})\right]_{x^{}=x}`$ (6)
$`=`$ $`{\displaystyle \frac{1}{\pi }}\underset{\epsilon 0^+}{lim}\mathrm{Im}{\displaystyle \underset{i=1}{\overset{N}{}}}\psi _i|\widehat{G}(x)\widehat{p}^n|\psi _i`$
and on its symmetrized form

$`S_n(x)`$ $`=`$ $`(i)^n\left[{\displaystyle \frac{^n}{r^n}}\rho ^{(l)}(x+r/2,xr/2)\right]_{r=0}`$ (9)
$`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle \frac{1}{2^n}}{\displaystyle \underset{\sigma =0}{\overset{n}{}}}\left(\begin{array}{c}n\\ \sigma \end{array}\right)\underset{\epsilon 0^+}{lim}\mathrm{Im}{\displaystyle \underset{i=1}{\overset{N}{}}}\psi _i|\widehat{p}^\sigma \widehat{G}(x)\widehat{p}^{n\sigma }|\psi _i.`$
Here we have set $`\mathrm{}=1`$ and $`m=1`$. For $`n=2`$ these definitions give rise respectively to twice the local kinetic energy $`T(x)`$ \[Eq. (6)\] and to the momentum flux density $`\mathrm{\Pi }(x)`$ \[Eq. (9)\]. Although the semiclassical limits and the integrals of these two functions coincide, in the quantum limit these are different functions: the kinetic energy density is invoked in density-functional theory, while the momentum flux density appears in the equations of generalized hydrodynamics.

The same method which was used in Sec. 2.1 to express the density profile in terms of the determinant of the inverse of the Green’s function operator can be extended to deal with the operators of the form $`\widehat{p}^\sigma \widehat{G}(x)\widehat{p}^{n\sigma }`$ . Again, for low values of $`n`$ one is reduced to evaluate with recursive methods the determinant of sparse matrices, which are always tridiagonal in the tails.

The kinetic energy and the momentum flux densities obtained with this procedure are shown in Fig. 2 for various number of fermions in 1D harmonic confinement. Oscillations along the profile and negative tails are found in the kinetic energy density, while a smooth positive profile is obtained for the momentum flux density (which however shows oscillations in its first derivative, see the inset of Fig. 2). Other properties of these curves are that the tails of the momentum flux density profile can be expressed analytically in terms of the von Weiszäcker surface energy $`(\mathrm{}^2/8m)|n(x)|^2/n(x)`$ , and that an approximate form of the kinetic energy density can be found using the exact density profile in the local-density expression $`T(x)=\pi ^2\mathrm{}^2[n(x)]^3/6m`$ .

### 2.3 Pair distribution function

Higher order correlations are described in the equal-times pair distribution function. We focus here on how the Pauli exclusion principle is reflected in the pair distribution function of a gas of $`N`$ noninteracting fermions in harmonic confinement. The same results apply to a 1D Bose gas with hard-core interactions, since the boson-fermion mapping holds for the pair distribution function .

The longitudinal pair distribution function of a quasi-1D Fermi gas is defined as

$$\rho _2^{(l)}(x_1,x_2)=n(x_1)n(x_2)F(x_1,x_2),$$
(10)
where the function $`F(x_1,x_2)`$ is given by

$$F(x_1,x_2)=\underset{i,j=1}{\overset{N}{}}\varphi _i^{}(x_1)\varphi _j(x_1)\varphi _j^{}(x_2)\varphi _i(x_2)$$
(11)
This function can be calculated by an extension of the Green’s function method used in the previous sections to evaluate particle, kinetic energy and momentum flux densities . We first rewrite Eq. (11) in terms of the Green’s function $`\widehat{G}(x)`$,

$$F(x_1,x_2)=\frac{1}{\pi ^2}\underset{\epsilon 0^+}{lim}\underset{i,j=1}{\overset{N}{}}\mathrm{Im}\psi _i|\widehat{G}(x_1)|\psi _j\mathrm{Im}\psi _j|\widehat{G}(x_2)|\psi _i;$$
(12)
we then use the property $`\mathrm{Im}A\mathrm{Im}B=(1/2)[\mathrm{Re}(AB^{})\mathrm{Re}(AB)]`$ to obtain the final expression

$$F(x_1,x_2)=\frac{1}{2\pi ^2}\underset{\epsilon 0^+}{lim}\left\{\mathrm{ReTr}\left[\widehat{G}_N(x_1)\widehat{G}_N^{}(x_2)\right]\mathrm{ReTr}\left[\widehat{G}_N(x_1)\widehat{G}_N(x_2)\right]\right\},$$
(13)
where $`\widehat{G}_N(x)`$ is the first $`N\times N`$ block of the matrix $`\widehat{G}(x)`$. In the case of harmonic confinement this can be evaluated again by making use of renormalization techniques and recursive methods: the evaluation of $`F(x_1,x_2)`$ is reduced to the calculation of the determinant of pentadiagonal $`N\times N`$ matrices, which can be factorized into the product of $`2\times 2`$ matrices.

With this Green’s function method we have calculated the longitudinal contribution to the equal-times pair distribution function for up to $`N=100`$ fermions without particular numerical efforts. Since this function presents in general a number of maxima of order $`N^2`$, for the sake of clarity we have shown in Fig. 3 the full result only for the case of $`N=4`$ fermions. We have plotted the pair distribution function in the center-of-mass $`R(x_1+x_2)/2`$ and relative $`r(x_1x_2)`$ coordinates: the effect of Pauli exclusion is observable in real space in the $`r`$ direction as a depression of $`\rho _2^{(l)}(x_1,x_2)`$ at short distances.

## 3 Dynamical properties

### 3.1 Collisionless regime

The dynamic structure factor for a non-interacting Fermi gas in the collisionless regime is given by

$$S(𝐤,\omega )=\underset{i,j}{}\left|d^3re^{i𝐤𝐫}\varphi _i^{}(𝐫)\varphi _j(𝐫)\right|^2f(\epsilon _i)[1f(\epsilon _j)]\mathrm{\hspace{0.17em}2}\pi \delta (\omega (\epsilon _j\epsilon _i)/\mathrm{}).$$
(14)
We focus here on quasi-1D excitations. The condition in this case is most stringent than for quasi-1D equilibrium properties (indeed the Paris experiment  satisfies the latter condition but not the former): the transverse excited states are not involved in the excitation processes when (i) the transverse component $`k_{}`$ of the momentum transfer vanishes, due to orthogonality of harmonic oscillator wave functions, or (ii) the energy transfer $`\omega `$ is smaller than the gap between the chemical potential and the first transverse excited state.

In the quasi-1D limit Eq. (14) reduces to a one-dimensional problem and the transverse-state wave function factorizes out. By evaluating the overlap integral in Fourier space and exploiting the properties of the Hermite polynomials , we obtain the expression for the dynamic structure factor as

$$𝒮(𝐤,h)=2\pi e^{k_{}^2/2\lambda }e^{k_x^2/2}\underset{i=\mathrm{max}\{Nh,0\}}{\overset{N1}{}}\frac{i!}{(i+h)!}\left(\frac{k_x^2}{2}\right)^h\left[L_i^h\left(k_x^2/2\right)\right]^2.$$
(15)
Here $`h`$ is an integer corresponding to a single-atom excitation of $`h`$ quanta of the harmonic oscillator, $`L_i^h(x)`$ is the $`i^{th}`$ generalized Laguerre polynomial of parameter $`h`$, $`\lambda `$ is the anisotropy parameter for the 3D harmonic oscillator and we have set $`\mathrm{}=1`$, $`m=1`$ and $`\omega _{ho}=1`$. The result for the spectrum is shown in Fig. 4 at two different values of the transferred wave vector $`k_x`$.

A good description for the spectrum is already given by the local-density approximation (LDA), defined as

$$S_{LDA}(k_x,\omega )=𝑑xn(x)S_{hom}(k_x,\omega ;\mu (x)).$$
(16)
The spatial dependence of the chemical potential $`\mu =k_F^2/2=\pi ^2n^2/2`$ is determined by the relation $`\mu (x)=\mu (n(x))`$ and $`n(x)=\sqrt{2Nx^2}/\pi `$; and the expression for the dynamic structure factor of the homogeneous system is $`S_{hom}(k_x,\omega )=\pi /k_xk_F`$ if $`|\omega _2(k_x)|<\omega <\omega _1(k_x)`$ and $`S_{hom}(k_x,\omega )=0`$ otherwise, with $`\omega _{1,2}=k_x^2/2\pm k_xk_F`$. The integral in Eq. (16) can be evaluated analytically to obtain

$$S_{LDA}(k_x,\omega )=\frac{2}{k_x}\left[\sqrt{2N\left(\frac{\omega }{k_x}\frac{k_x}{2}\right)^2}\sqrt{2N\left(\frac{\omega }{k_x}+\frac{k_x}{2}\right)^2}\right].$$
(17)
From Fig. 4 we see that the LDA description captures the main features of the spectrum. The presence of the external harmonic potential significantly modifies the spectrum with respect that of the homogeneous system. Due to the boson-fermion mapping  the spectrum represented in Fig. 4 is also valid for the Tonks gas.

### 3.2 Collisional regime

We turn now to the long-wavelength limit and we investigate the analogue of sound-wave propagation. Due to the presence of the external confinement, the collective excitations are quantized. We evaluate here the spectrum and the expression for the density fluctuations in the linear regime.

The equation of motion for the density profile, as can be derived exactly from the equations of generalized hydrodynamics , takes the form

$$m_t^2n(x,t)=_x^2\mathrm{\Pi }(x,t)+_x\left[n(x,t)_xV(x,t)\right].$$
(18)
We define the collisional regime by requiring that the momentum flux density $`\mathrm{\Pi }(x,t)`$ depends locally on the particle density and the velocity field through the relation

$$\mathrm{\Pi }(x,t)=\mathrm{}^2\pi ^2n^3(x,t)/3m+mn(x,t)v^2(x,t)/2.$$
(19)
In the linear regime equations (18) and (19) lead to a closed equation for the density fluctuation $`\delta n(x,t)`$. In the homogeneous limit we obtain a phonon dispersion relation with velocity $`c=\mathrm{}k_F/m`$.

In the case of 1D harmonic confinement in order to find a discrete spectrum we then need to assume a generalized form of the Thomas-Fermi approximation (TFA), which reads

$$\delta \mathrm{\Pi }(x,t)2\frac{\delta t}{\delta n}\delta n(x,t)=2[\mu V_{ext}(x)]\delta n(x,t)$$
(20)
where $`t(n)`$ is the kinetic energy density and the second equality follows from the Euler equation for density functional theory. We assume that Eq. (20) holds also outside the classical radius $`X_F=\sqrt{2m\mu /\omega _{ho}^2}`$: this hypothesis is necessary since the solutions of the linearized Eqs. (18) and (19) diverge at the classical boundary as $`|x^2X_F^2|^{1/2}`$ and thus we cannot obtain the dispersion relation by imposing that the solution should vanish at $`x=X_F`$, as in the 3D case . By matching the solutions $`\delta n(x)\sqrt{|x^2X_F^2|}`$ inside and outside the classical radius we obtain the spectrum

$$\omega =n\omega _{ho},$$
(21)
with corresponding solutions inside the classical radius

$$\delta n_{in}(z)=\frac{\mathrm{cos}(n\mathrm{arccos}(z))}{\sqrt{1z^2}}$$
(22)
and outside the classical radius

$$\delta n_{out}(z)=\frac{\left(|z|\sqrt{z^21}\right)^n}{\sqrt{z^21}}.$$
(23)
Here we have set $`z=x/X_F`$. The divergence of the solutions at $`x=X_F`$ is unphysical, and is a consequence of our approximations: the linearized TFA solutions are not valid around the classical turning points.

### 3.3 Non-Linear Schrödinger equation

A possibility to go beyond the TFA is to describe the collisional Fermi gas through the non-linear Schrödinger equation

$$i\mathrm{}_t\mathrm{\Phi }(x,t)=\left[\frac{\mathrm{}^2}{2m}_x^2+V(x,t)+\frac{\pi ^2\mathrm{}^2}{2m}|\mathrm{\Phi }(x,t)|^4\right]\mathrm{\Phi }(x,t).$$
(24)
Here, $`\mathrm{\Phi }(x,t)`$ is normalized to the number $`N`$ of particles in the trap. According to Kolomeisky et al. , Eq. (24) describes the motion of a 1D gas of impenetrable bosons in a long-wavelength approach. Although Eq. (24) has been shown to fail in the description of the dynamics of the phase of the Bose gas , in the linear regime this equation describes the collisional dynamics of the Fermi gas. Indeed, by setting $`\mathrm{\Phi }(x,t)=\sqrt{n(x,t)}\mathrm{exp}[i\varphi (x,t)]`$ one can transform Eq. (24) into hydrodynamic equations for the density $`n(x,t)`$ and the phase $`\varphi (x,t)`$ of the fluid; linearization and neglect of kinetic energy term $`_x^2\sqrt{n(x,t)}/2m\sqrt{n(x,t)}`$ leads to an expression which coincides with the linear form of Eqs. (18) and (19).

In the case of a harmonic external potential a numerical solution for the low-lying collective modes of Eq. (24) has been obtained elsewhere. The numerical results are in agreement with the analytical spectrum (21) and with the main features of the density fluctuations profiles \[Eqs. (22) and (23)\] .

It is striking that the solution of Eq. (24) for the collective modes in harmonic confinement always yields the free harmonic-oscillator spectrum (21), independently of the strength of the interactions. In the case of the quadrupole mode, this result can be related to an underlying scaling symmetry of the Hamiltonian giving rise to Eq. (24). This is the one-dimensional analogue  of the hidden symmetry already noticed by Pitaevskii and Rosch for the 2D Gross-Pitaevskii equation : in the Hamiltonian the interaction-potential term

$$v(xx^{})=\frac{\mathrm{}^2\pi ^2}{3m}n(x)\delta (xx^{})$$
(25)
scales under space dilatations $`x\lambda x`$ in the same fashion as the kinetic-energy term, and this property in the case of external harmonic confinement leads to a closed equation of motion for the expectation value of the quadrupole-mode operator $`I=m\widehat{x}^2/2`$, expressed only in terms of the total energy $`E`$:

$$_t^2I=4\omega _{ho}^2I+4E.$$
(26)
This yields $`\omega =2\omega _{ho}`$ for the frequency of the quadrupole mode. It is interesting to compare this value with what is obtained from the 1D Gross-Pitaevskii equation , namely $`\omega =\sqrt{5/2}\omega _{ho}`$.

## 4 Summary and concluding remarks

In summary, we have studied the equilibrium and dynamical properties of a one-dimensional Fermi gas under harmonic confinement. We have developed a method which allows to evaluate the exact equilibrium profiles of the system for up to a large number of fermions. The density profiles show a shell structure which is very prominent in this reduced dimensionality. We have calculated the excitation spectrum of the Fermi gas both in the collisionless regime, where an analytic exact expression can be obtained, and in the collisional regime, where we have given a solution in the generalized Thomas-Fermi approximation. The spectra obtained in the two regimes coincide – this is the analogue of the coincidence of the velocities of zero and first sound in the 1D homogeneous Fermi gas. The results here presented for the density profiles, pair distribution function and excitation spectrum describe as well the 1D Bose gas in the Tonks regime. The measurement of excitation frequency of the quadrupole modes should be a sensible probe of the Tonks-gas phase.

We thank Professor Tosi for useful discussions and encouragement. This research is supported by INFM under the project PRA2001.

## References


format: text

text:
# I Introduction

## I Introduction

While the good amount of evidences for the neutrino mass and the lepton flavor mixing have been accumulated , we still lack observational indications of how large is the absolute mass of the neutrinos. To our understanding to date it may show up in only a few places, the single  or the double beta decay  experiments as well as future cosmological observations . Other potential possibilities for hints of absolute mass of neutrinos include Z-burst interpretation of highest energy cosmic rays .

Among these various experimental possibilities the neutrinoless double beta decay experiments seems to have relatively higher sensitivities. The most stringent bound on effective mass parameter $`m_{\beta \beta }`$ (see eq. (3) for definition) is now $`m_{\beta \beta }<0.35`$ eV, which comes from Heidelberg-Moscow group . Furthermore, a wide variety of proposals for future facilities as well as ongoing attempt with greater sensitivies are actively discussed. They include NEMO , GENIUS , CUORE , MOON , XMASS , and EXO  projects. These high-sensitivety experiments open the enlighting possibility of discovering neutrinoless double beta decay events, not just placing an upper bound on $`m_{\beta \beta }`$ by its nonobservation. Therefore, it is of great importance to completely understand what kind of informations can be extracted if such discovery is made.

We discuss in this paper in a generic three flavor mixing framework the constraints on neutrino masses and mixing by positive observation (as well as nonobservation) of neutrinoless double beta decay. The constraints imposed on neutrino mixing parameters by neutrinoless double beta decay have been discussed by many authors. They include the ones in early epoch , those in ”modern era” in which real constraints on solar mixing parameters are started to be extracted , and the ones in ”post-modern era” where the analyses are performed in a complehensive manner in the framework of generic three flavor neutrino mixing .

In a previous paper, we have made a final step in the series of analyses by proposing a way of expressing the constraints solely in terms of observables in single and double beta decay . By using the framework, we discussed the possibility of placing lower bound on $`|U_{e3}|^2`$ assuming positive observation in direct mass measurement in single beta decay and an upper limit on $`m_{\beta \beta }`$ in double beta decay experiments. It is a natural and logical step for us to examine next the alternative case of positive observation of neutrinoless double beta decay events.

Timely enough, an evidence of the neutrinoless double beta decay has just been reported by Klapdor-Kleingrothaus and collaborators . Since the confidence levels of the claimed evidence are about 2 and 3 $`\sigma `$ in Bayesian and Partcle Data Group methods, respectively, we must wait for confirmation by further data taking, or by other groups to conclude that neutrinos are Majorana particles. Nevertheless, we feel that the peak in the relevant kinematic region in their experiments is too prominent to be simply ignored.

As will become clear as we proceed it is essential to combine the constraint on $`|U_{e3}|^2=s_{13}^2`$ imposed by the reactor experiments . One of the key points in our subsequent discussion is that the double beta and the reactor bounds cooperate to produce a stringent constraint on absolute mass scale of neutrinos and the mixing angle $`\theta _{12}`$ which is responsible for the solar neutrino problem.

## II Constraints from neutrinoless double beta decay

Let us start by defining our notations. We use throughout this paper the standard notation of the MNS matrix :

$$U=\left[\begin{array}{ccc}c_{12}c_{13}& s_{12}c_{13}& s_{13}e^{i\delta }\\ s_{12}c_{23}c_{12}s_{23}s_{13}e^{i\delta }& c_{12}c_{23}s_{12}s_{23}s_{13}e^{i\delta }& s_{23}c_{13}\\ s_{12}s_{23}c_{12}c_{23}s_{13}e^{i\delta }& c_{12}s_{23}s_{12}c_{23}s_{13}e^{i\delta }& c_{23}c_{13}\end{array}\right].$$
(1)
Using the notation, the observable in neutrinoless double beta decay experiments can be expressed as

$`m_{\beta \beta }`$ $`=`$ $`\left|{\displaystyle \underset{i=1}{\overset{3}{}}}m_iU_{ei}^2\right|`$ (2)
$`=`$ $`\left|m_1c_{12}^2c_{13}^2e^{i\beta }+m_2s_{12}^2c_{13}^2e^{+i\beta }+m_3s_{13}^2e^{i(3\gamma 2\delta )}\right|,`$ (3)
where $`m_i`$ (i=1, 2, 3) denote neutrino mass eigenvalues, $`U_{ei}`$ are the elements in the first low of the MNS matrix, and $`\beta `$ and $`\gamma `$ are the extra CP-violating phases characteristic to Majorana neutrinos , for which we use the convention of Ref. .

We define the neutrino mass-squared difference as $`\mathrm{\Delta }m_{ij}^2m_j^2m_i^2`$ in this paper. In the following analysis, we must distinguish the two different neutrino mass patterns, the normal ($`\mathrm{\Delta }m_{23}^2>0`$) vs. inverted ($`\mathrm{\Delta }m_{23}^2<0`$) mass hierarchies. We use the convention that $`m_3`$ is the largest (smallest) mass in the normal (inverted) mass hierarchy so that the angles $`\theta _{12}`$ and $`\theta _{23}`$ are always responsible for the solar and the atmospheric neutrino oscillations, respectively. We therefore often use the notations $`|\mathrm{\Delta }m_{23}^2|\mathrm{\Delta }m_{atm}^2`$ and $`\mathrm{\Delta }m_{12}^2\mathrm{\Delta }m_{}^2`$ to emphasize that they are experimentally measurable quantities. Because of the hierarchy of mass scales, $`\mathrm{\Delta }m_{}^2/\mathrm{\Delta }m_{atm}^21`$, $`\mathrm{\Delta }m_{12}^2`$ can be made always positive as far as $`\theta _{12}`$ is taken in its full range \[0, $`\pi /2`$\].

In order to derive constraint on mixing parameters we need the classification.

$$\text{Case I:}\text{ }\left|m_1c_{12}^2c_{13}^2e^{i\beta }+m_2s_{12}^2c_{13}^2e^{+i\beta }\right|m_3s_{13}^2$$
(4)
$$\text{Case II:}\text{ }\left|m_1c_{12}^2c_{13}^2e^{i\beta }+m_2s_{12}^2c_{13}^2e^{+i\beta }\right|m_3s_{13}^2$$
(5)
However, examination of the case II reveals that it does not lead to useful bounds. Therefore, we only discuss the case I in the rest of this paper.

### A Joint constraint by upper bounds on $`m_{\beta \beta }`$ and reactor experiments

Since we try to utilize the experimental upper bound on $`m_{\beta \beta }`$, $`m_{\beta \beta }m_{\beta \beta }^{max}`$, we derive the lower bound on $`m_{\beta \beta }`$. It can be obtained in the following way;

$`m_{\beta \beta }`$ $``$ $`c_{13}^2\left|(m_1c_{12}^2+m_2s_{12}^2)\mathrm{cos}\beta i(m_1c_{12}^2m_2s_{12}^2)\mathrm{sin}\beta \right|m_3s_{13}^2`$ (6)
$`=`$ $`c_{13}^2\sqrt{m_1^2c_{12}^4+m_2^2s_{12}^4+2m_1m_2c_{12}^2s_{12}^2\mathrm{cos}2\beta }m_3s_{13}^2.`$ (7)
Noticing that the right-hand-side (RHS) of (7) has a minimum at $`\mathrm{cos}2\beta =1`$, we obtain the inequality

$$m_{\beta \beta }c_{13}^2\left|m_1c_{12}^2m_2s_{12}^2\right|m_3s_{13}^2.$$
(8)
We note that the RHS of (8) is a decreasing function of $`s_{13}^2`$, and hence takes a minimum value for the maximum value of $`s_{13}^2`$ which is allowed by the limit placed by the reactor experiments . We denote the maximum value as $`s_{13}^2(\text{CH})`$ throughout this paper. Numerically, $`s_{13}^2(\text{CH})0.03`$. (While the precise value of the CHOOZ constraint actually depends upon the value of $`\mathrm{\Delta }m_{atm}^2`$ , we do not elaborate this point in this paper.) Using the constraint we obtain

$$m_{\beta \beta }^{max}m_{\beta \beta }\left|m_1c_{12}^2m_2s_{12}^2\right|\left(m_3+\left|m_1c_{12}^2m_2s_{12}^2\right|\right)s_{13}^2(\text{CH}).$$
(9)
It can be rewritten as the bound on $`\mathrm{cos}2\theta _{12}=\mathrm{cos}2\theta _{}`$ as

$$\frac{m_2m_1}{m_2+m_1}\frac{m_{\beta \beta }^{max}+m_3s_{13}^2(\text{CH})}{\frac{1}{2}(m_2+m_1)c_{13}^2(\text{CH})}\mathrm{cos}2\theta _{12}\frac{m_2m_1}{m_2+m_1}+\frac{m_{\beta \beta }^{max}+m_3s_{13}^2(\text{CH})}{\frac{1}{2}(m_2+m_1)c_{13}^2(\text{CH})},$$
(10)
where $`c_{13}^2(\text{CH})1s_{13}^2(\text{CH})`$.

### B Joint constraint by lower bounds on $`m_{\beta \beta }`$ and reactor experiments

A positive observation of neutrinoless double beta decay will lead to the experimental lower bound on $`m_{\beta \beta }`$, $`m_{\beta \beta }m_{\beta \beta }^{min}`$, which we use to place new bound on neutrino mixing parameters. Toward the goal we note, similarly as (7), that

$$m_{\beta \beta }c_{13}^2\sqrt{m_1^2c_{12}^4+m_2^2s_{12}^4+2m_1m_2c_{12}^2s_{12}^2\mathrm{cos}2\beta }+m_3s_{13}^2,$$
(11)
whose RHS is maximized by taking $`\mathrm{cos}2\beta =+1`$ and $`s_{13}^2=s_{13}^2(\text{CH})`$ in the last term and $`c_{13}^2=1`$ in front of the square root. (A more refined treatment entails the same excluded region.) One can then derive an inequality similar to (9);

$$m_{\beta \beta }^{min}m_{\beta \beta }\left(m_1c_{12}^2+m_2s_{12}^2\right)+m_3s_{13}^2(\text{CH}).$$
(12)
By rewriting (12) we obtain the other upper bound on $`\mathrm{cos}2\theta _{12}`$;

$$\mathrm{cos}2\theta _{12}\frac{m_2+m_1}{m_2m_1}\frac{m_{\beta \beta }^{min}m_3s_{13}^2(\text{CH})}{\frac{1}{2}(m_2m_1)}.$$
(13)
To summarize, we have derived in this section the two kinds of upper bound on $`\mathrm{cos}2\theta _{12}`$ (lower bound for $`\mathrm{cos}2\theta _{12}<0`$) by using the assumed experimental constraint $`m_{\beta \beta }^{min}m_{\beta \beta }m_{\beta \beta }^{max}`$ imposed by neutrinoless double beta decay experiments.

## III Constraints expressed by experimental observables

We rewrite the bounds on solar mixing angle in terms of measurable quantities. Toward the goal we note that three neutrino masses $`m_i`$ (i=1,2,3) can be expressed by the two $`\mathrm{\Delta }m^2`$ and a remaining over-all scale $`m_H`$. We assign $`m_H`$ to the mass of the highest-mass state, $`m_3`$ in the normal mass hierarchy ($`\mathrm{\Delta }m_{23}^2>0`$), and $`m_2`$ in the inverted mass hierarchy ($`\mathrm{\Delta }m_{23}^2<0`$), respectively.

We have argued in our previous paper  that in a reasonable approximation one can regard $`m_H`$ as the observable $`m_\beta `$ in direct mass measurements in single beta decay experiments.***While we used the linear formula derived by Farzan, Peres and Smirnov
$$m_\beta =\frac{_{j=1}^nm_j|U_{ej}|^2}{_{j=1}^n|U_{ej}|^2}$$
(14)  with $`n`$ being the dimension of the subspace of (approximately) degenerate mass neutrinos, this point remains valid even if we use an alternative quadratic expression . We have noticed that the identification is exact in two extreme cases of degenerate and hierarchical mass spectra. Then, the three mass eigenvalues of neutrinos can be represented solely by observables; $`\mathrm{\Delta }m_{atm}^2`$, $`\mathrm{\Delta }m_{}^2`$, and $`m_\beta `$ in a good approximation.

In each neutrino mass pattern, we have the expressions of three mass eigenvalues:

Normal mass hierarchy ($`\mathrm{\Delta }m_{23}^2>0`$);

$$m_1=\sqrt{m_H^2\mathrm{\Delta }m_{atm}^2\mathrm{\Delta }m_{}^2},m_2=\sqrt{m_H^2\mathrm{\Delta }m_{atm}^2},m_3=m_Hm_\beta .$$
(15)
Inverted mass hierarchy ($`\mathrm{\Delta }m_{23}^2<0`$);

$$m_1=\sqrt{m_H^2\mathrm{\Delta }m_{}^2},m_2=m_Hm_\beta ,m_3=\sqrt{m_H^2\mathrm{\Delta }m_{atm}^2}.$$
(16)
It is instructive to work out the form of constraint in the degenerate mass approximation, $`m_i^2m^2\mathrm{\Delta }m_{atm}^2,\mathrm{\Delta }m_{}^2`$. It is easy to show that in the degenerate mass limit the bound (10) becomes

$$|\mathrm{cos}2\theta _{12}|\mathrm{sec}^2\theta _{13}(\text{CH})\left[\frac{m_{\beta \beta }^{max}}{m_\beta }+s_{13}^2(\text{CH})\right].$$
(17)
On the other hand, the bound (11) gives the inequality $`m_\beta m_{\beta \beta }^{min}`$ in the degenerate mass limit. (To show this one may go back to (11), rather than using (13).)

## IV Analysis of the Double Beta-Reactor Joint Constraints

We analyze in this section the joint constraints derived in the foregoing sections and try to extract the implications. Let us start by examining the case of recent obsevation announced in  which gives rise to 0.11 eV $`m_{\beta \beta }`$ 0.56 eV and 0.05 eV $`m_{\beta \beta }`$ 0.84 eV if $`\pm 50`$ % uncertainty of the nuclear matrix elements are considered, each at 95 % CL. In Fig. 1 we present on $`m_\beta `$ \- $`\mathrm{cos}2\theta _{12}`$ plane the constraint (10) by the thick solid lines (solid line) and (13) by the thick dashed line (dashed line) for cases with (without) uncertainty of the nuclear matrix elements, respectively. The regions surrounded by these lines are allowed. The slope of $`m_\beta `$-dependence of (13) is so large that the dashed line looks like a vertical line, which implies the inequality $`m_\beta m_{\beta \beta }^{min}`$. We have derived it earlier in the degenerate mass limit, but it is generically true if $`\mathrm{\Delta }m_{}^2`$ is smaller than other relevant mass squared scales. Only the case of normal mass hierarchy ($`\mathrm{\Delta }m_{23}^2>0`$) is shown in Fig. 1; the case of inverted hierarchy ($`\mathrm{\Delta }m_{23}^2<0`$) gives an almost identical result except for a slight shift of the dashed line toward smaller $`m_\beta `$ by $``$ 10 %.

Superimposed in Fig. 1 are the 95 % CL allowed regions of $`\mathrm{cos}2\theta _{12}`$ for the large mixing angle (LMA) MSW solution (indicated by the shaded region between thin solid lines) and the low (LOW) MSW solution (indicated by the shaded region between thin dashed lines) of the solar neutrino problem . There are several up to date global analyses of the solar neutrino data with similar results of allowed region of mixing parameters . Therefore, we just quote the result obtaind by Krastev and Smirnov in the last reference in .

Figure 1 illustrates that for a given value of $`\mathrm{cos}2\theta _{12}`$ the single beta decay observable $`m_\beta `$ has to fall into a region bounded by $`m_\beta ^{min}m_{\beta \beta }^{min}`$ and $`m_\beta ^{max}`$, which are dictated by (13) and (10), respectively. Thus, we have a clear prediction for direct mass measurements using a single beta decay with observation of double beta decay events. With use of the numbers given in , for example, the observable $`m_\beta `$ must fall into the region 0.05 eV $`m_\beta `$ 2 eV (0.11 eV $`m_\beta `$ 1.3 eV) with (without) uncertainty of nuclear matrix elements at the best fit parameters of the LMA MSW solution. (The best fit value is $`\mathrm{tan}^2\theta _{12}=0.35`$, or $`\mathrm{cos}2\theta _{12}=0.48`$ in the last reference in .) Within the allowed region the cancellation between three mass eigenstates can take place for appropriate values of Majorana phases that allow (typically) a factor of 2-3 larger values of $`m_\beta `$ compared with the measured value of $`m_{\beta \beta }`$. At around maximal mixing ($`\mathrm{cos}2\theta _{12}0`$), which is allowed by 95 % CL in the LOW solution, the cancellation is so efficient that much larger values of $`m_\beta `$ is allowed. Therefore, there are still ample room for hot dark matter mass neutrinos both in the LMA and the LOW solutions.

It should be emphasized that a finite value of $`m_{\beta \beta }`$ does imply a lower bound on $`m_\beta `$, as indicated in Fig. 1; a vanishingly small $`m_\beta `$ cannot be consistent with finite $`m_{\beta \beta }`$ in double beta decay experiments. The sensitivity of the proposed KATRIN experiment is expected to extend to $`m_\beta 0.3`$ eV . On the other hand, the present 68 % CL limit quoted in  without nuclear element uncertainty is 0.28 eV $`m_{\beta \beta }`$ 0.49 eV. Therefore, if the limit is further tightened by additional data taking in the future, both experiments can become inconsistent, giving an another opportunity of cross checking.

In Fig. 2, we demonstrate the approximate scaling relation obeyed by the constraint (10) by taking $`m_\beta /m_{\beta \beta }^{max}`$ as the abscissa in a wide range of the $`m_{\beta \beta }^{max}`$ in degerarate mass region, 0.1 eV $`<m_{\beta \beta }^{max}<`$ 1 eV. The scaling is exact in the degenerate mass limit as shown in (17). The relation is useful to estimate the allowed region of $`m_\beta `$ for a given value of $`m_{\beta \beta }^{max}`$ which is not explicitly examined in this paper.

The most stringent bound to date on $`m_\beta `$ is from the Mainz collaboration , $`m_\beta 2.2`$ eV (95 % CL). (A similar bound $`m_\beta 2.5`$ eV (95 % CL) is derived by the Troitsk group .) As we can see in Fig. 2 that the double beta bound with the CUORE sensitivity region $`m_{\beta \beta }<0.3`$ eV  becomes stronger than the Mainz bound for the LMA MSW solution but not for the LOW MSW solution in their 95 % CL regions.

In Fig. 3, we present the similar allowed regions for hypothetical discovery of neutrinoless double beta decay events which would produce the experimental bounds 0.01 eV $`m_{\beta \beta }`$ 0.03 eV. It is to examine how the constraint changes in some other situation of discovery with different mass parameter ranges. We note that even such deep region of sensitivity will be explored by several experiments .

For this case, the bounds for the normal and the inverted mass hierarchies start to split as shown in Fig. 3. In the case of inverted mass hierarchy the lower bound on $`m_\beta `$ is replaced by the trivial bound $`m_\beta \sqrt{\mathrm{\Delta }m_{atm}^2}`$ which is more restrictive. The latter is indicated by the dash-dotted line in Fig. 3b. It is also evident that the constraint from double beta decay is so stringent that the limit on $`m_\beta `$ is tightened to be $`m_\beta <`$ 0.2 eV for the LMA MSW solution.

In conclusion, we have demonstrated in this and the previous papers the mutual intimate relationship between observation and/or nonobservation in single beta decay and neutrinoless double beta decay experiments. We hope that it stimulates even richer future prospects not only in double beta decay experiments but also in direct mass measurements using single beta decay.

Finally, some remarks are in order:

(1) If the LMA MSW solution is the case and if the KamLAND experiment  that just started data taking can measure $`\mathrm{cos}2\theta _{12}`$ within 10 % level accuracy, the upper limit of $`m_\beta `$ can be accurately determined with $``$ 20 % accuracy.

(2) In this paper we have derived constraints imposed on neutrino mixing parameters by observation of neutrinoless double beta decay events and the CHOOZ reactor bound on $`|U_{e3}|^2`$ in the generic three flavor mixing framework of neutrinos. Suppose that neutrinoless double beta decay events are confirmed to exist and the single beta decay experiments detect neutrino mass outside the region of the bound derived in this paper. What does it mean? It means either that double beta decay would be mediated by some mechanisms different from the usual one with Majorana neutrinos, or the three flavor mixing framework used in this paper is too tight.

Nore added:

After submitting the first version of our paper to the electronic archive, we became aware of the works which address relatively model-independent implication of the results reported in , or critically comment on the interpretation of the events. References  and  are the incomplete lists of them.

###### Acknowledgements.

We thank Professor Klapdor-Kleingrothaus for calling our attention to the reference . This work was supported by the Grant-in-Aid for Scientific Research in Priority Areas No. 12047222, Japan Ministry of Education, Culture, Sports, Science, and Technology.


text:
# Density dependent spin polarisation in ultra low-disorder quantum wires

## Abstract

There is controversy as to whether a one-dimensional (1D) electron gas can spin polarise in the absence of a magnetic field. Together with a simple model, we present conductance measurements on ultra low-disorder quantum wires supportive of a spin polarisation at $`B=0`$. A spin energy gap is indicated by the presence of a feature in the range $`0.50.7\times 2e^2/h`$ in conductance data. Importantly, it appears that the spin gap is not static but a function of the electron density. Data obtained using a bias spectroscopy technique are consistent with the spin gap widening further as the Fermi-level is increased.

In the presence of strong exchange coupling, electrons can spin polarise in the absence of an applied magnetic field. Such a scenario is predicted for a variety of different systems including one dimensional (1D) ballistic quantum wires wang1 ; Goldgpar , the two dimensional (2D) electron gas 2DEGpol , three dimensional (3D) metal nanowires zabalaPRL  and circular quantum dots koskinenPRL97 . In the case of 1D, interactions become increasingly important at low densities and models such as the Tomanaga-Luttinger liquid theory Luttinger  are required to describe them. Despite the large exchange energy present in low-density 1D systems there are strict theoretical arguments against magnetic ordering leibmattis  and the notion of a 1D spin-polarised ground state remains the subject of wide debate, in particular since the important experimental results of Thomas et al. Thomasspin1st . Here we present experimental results taken on semiconductor quantum wires in zero magnetic field that provide strong evidence in favor of a spin energy gap developing in the 1D region. This density-dependent energy gap between spin-up and spin-down electrons is revealed in transport measurements as an anomalous conductance feature in the range $`0.50.7\times 2e^2/h`$.

A feature near $`0.7\times 2e^2/h`$ can be seen in some of the earliest 1D ballistic transport measurements on quantum point contacts Wharam ; vanWeesQPC . In 1996 Thomas et al. Thomasspin1st  revealed that the anomalous feature was related to spin by showing that it evolves smoothly into the Zeeman spin-split level at $`0.5\times 2e^2/h`$ with an in-plane magnetic field. Since that time experimental studies have concentrated on the behavior of the anomaly as a function of temperature, source-drain bias, magnetic field, thermopower, wire length and density Thomasint ; BKQWAPL ; ct-liang ; kristensenPRB ; newthomas ; Pyshkin ; appleyard ; reilly . Together with these investigations, numerous mechanisms to explain the origin of the conductance feature have been proposed schmeltzer ; flambaum ; spivak ; wingreen1 ; sushkov ; bruus2 ; wingreen2 . Amongst the most compelling of these models is the notion of Fermi-level pinning in the presence of a static spin energy gap bruus2 ; wingreen2 .

In this work we propose, and present supportive experimental data for, an alternative simple phenomological model that appears to explain the characteristic details of the $`0.7`$ feature by means of a density dependent spin gap arising in the region of the quantum wire. Key differences exist between this simple model and other explanations based on Fermi level pinning. As depicted in Fig.(1), the spin levels are degenerate when the sub-band is empty, with the spin gap opening up as the carrier density in the 1D sub-band is increased. Consequently this model also explains the absence of conductance plateaus at $`0.25\times 2e^2/h`$ in the presence of a finite source - drain (S-D) bias.

In contrast to the case of a static spin polarisation, the sub-bands remain spin degenerate until they are populated, after which the spin gap opens with increasing density. This behavior, suggestive of many-body interactions, is consistent with a spin polarisation driven by exchange as predicted by Wang and Berggren wang2 . In their calculations the polarisation weakens as higher sub-bands are populated (see Fig.1). The non-linear dependence of Fermi energy $`E_F`$ on density or gate voltage in Fig.1 is a consequence of the singularity in the 1D density of states, $`\rho E^{1/2}`$ wingreen1 .

Our proposed model is consistent with conductance data we have obtained on GaAs/AlGaAs quantum wires free from the disorder associated with modulation doping BKQWAPL . Although we illustrate our simple model via comparison with this data, the model is not limited to these samples but appears to be consistent with the key published results Thomasspin1st ; Wharam ; vanWeesQPC ; Thomasint ; BKQWAPL ; ct-liang ; newthomas ; kristensenPRB ; Pyshkin ; appleyard ; reilly .

Turning now to the experimental data, we note that the presence and shape of conductance anomalies depend on the rate at which the spin gap opens with 1D density and are likely to be sample dependent. Figure 2 depicts the three main scenarios that may arise. Scenario (I) occurs if spin splitting takes place quickly with increasing density, so that an appreciable energy gap develops in comparison to the thermal energy $`kT`$. In this case we see a fully resolved spin-split plateau near $`0.5\times 2e^2/h`$ in linear response conductance $`G`$ and no feature near $`0.7\times 2e^2/h`$. This is shown on the left of Figure 2(I), for a quantum wire of length $`l=0.5\mu m`$ at $`T`$ = 50mK. The right side of the figure shows the dependence of the differential conductance ($`di/dv`$) with finite source-drain bias, where the thick lines represent conductance plateaus. Due to an averaging of the conductance at the chemical potential of the source $`S`$ and drain $`D`$, half-plateaus at 0.5 and 1.5 $`\times 2e^2/h`$ occur at finite bias when the two potentials differ by one sub-band patel ; moreno .

In a large applied magnetic field the Zeeman energy lifts the spin degeneracy, such that plateaus at $`0.5,1.5,2.5\mathrm{}\times 2e^2/h`$ are resolved in the the conductance Wharam . The simultaneous application of a finite S - D bias then produces additional quarter plateaus at $`0.25,0.75,1.25,1.75\mathrm{}\times 2e^2/h`$ patelparrab . However, in the case of a density dependent spin-gap at $`B`$ = 0, the plateaus at $`0.25`$ and $`1.25\times 2e^2/h`$ will be absent, since the 1D electron density is not large enough to appreciably open the spin gap.

Scenario (II) considers the case where the thermal energy $`kT`$ is comparable to the spin gap. In this case no feature near $`0.5\times 2e^2/h`$ will be resolved as $`E_F`$ crosses the band edges. Instead, as $`E_F`$ approaches the upper spin band-edge, the spin-gap continues to open so that the number of electrons which thermally populate the upper spin band remains approximately constant. A quasi-plateau near $`0.7\times 2e^2/h`$ therefore occurs due to the simultaneous increase of both $`E_F`$ and the upper spin-split sub-band, (shown in the data for the same $`l=0.5\mu m`$ wire at $`T`$ = 4.2K). In this model the quasi-plateau can occur anywhere in the range $`0.51.0\times 2e^2/h`$, as has been observed experimentally newthomas ; reilly .

The right side of Fig. 2(II) illustrates the behavior of the differential conductance at elevated temperatures. In contrast to the low temperature case of scenario (I), the feature remains close to $`0.75\times 2e^2/h`$ even when the S - D bias is close to zero (Fig. 2(II) d & e).

Scenario (III) illustrates the case for where the spin-splitting is weak and only grows slowly with increasing density. At low temperatures a feature near $`0.5\times 2e^2/h`$ is absent if the spin gap remains small in comparison to $`kT`$. We illustrate this scenario with data taken on a $`l`$ = 0 quantum wire at $`T`$ = 50mK. Although there is no evidence for a spin gap at zero S - D bias, the spin gap can still be observed as a feature at $`0.75\times 2e^2/h`$ in the differential conductance as shown in the schematic on the right. This is because the splitting is too small to be resolved at low densities where $`G<2e^2/h`$ at $`V_{SD}`$ = 0. However, a moderate S - D bias evolves the $`1.0\times 2e^2/h`$ plateau into a feature near $`0.75\times 2e^2/h`$ as the source and drain potentials differ by one spin sub-band (Fig. 2(II)e). Such behavior has been observed by Kristensen et al. kristensenPRB ; foot2 .

We now present conductance data taken on two different samples which support our model. The devices are fabricated from ultra low-disorder GaAs/AlGaAs heterostructures with electron mobilities in the range $`46\times 10^6cm^2V^1s^1`$. Using electron beam lithography the top layer of the heterostructure is sectioned into three separately controllable gates with the middle or top gate being biased positively ($`V_T`$) to control both the density in the 2D reservoirs ($`n_{2D}=1.57.5\times 10^{11}`$cm<sup>-2</sup>) and in the 1D channel BKQWAPL . The side gates are negatively biased ($`V_S`$) and simultaneously control the 1D density and the transverse potential. Careful control of all three independent gates allows the $`0.7\times 2e^2/h`$ feature to be studied as a function of both density and potential profile. We find that the feature tends towards $`0.5\times 2e^2/h`$ with increasing top gate bias and length of the 1D region reilly .

Figure 3(a) shows data taken on a quantum wire of length $`l=0.5\mu m`$ at $`T`$ = 4.2K and $`T`$ = 50mK. The data taken at $`T`$ = 50mK show a characteristic evolution towards fully resolved spin-splitting with the $`0.7`$ feature moving closer to $`0.5\times 2e^2/h`$ with increasing top gate bias (right to left). In the context of the model this evolution is consistent with the spin gap opening more rapidly with 1D density $`n`$, for larger $`V_T`$, (i.e. moving from scenario III to scenario I with increasing $`V_T`$) although the exact mechanism behind the dependence on $`V_T`$ is presently unknown foot . In particular, the left-most traces from Fig. 3(a) are consistent with scenario (I) at $`T`$ = 50mK and scenario (II) at $`T`$ = 4.2K, where the feature has risen from $`0.55`$ to $`0.7\times 2e^2/h`$.

At $`T`$ = 4.2K the position of the feature does not evolve with top gate bias (2D density) but remains close to $`0.65\times 2e^2/h`$ at these higher temperatures, as both of the 1D spin bands are populated (c.f. Fig 2.(II)), and the position of the feature is rather insensitive to small changes in the spin-gap.

The shape of the feature will depend on the slope of the Fermi function as it crosses the second spin-band edge: at high temperatures the broad Fermi function produces a broad quasi-plateau. At low temperatures, depending on the size of the spin gap, the sharp Fermi function will either produce a small but sharp feature at $`0.5\times 2e^2/h`$ (strong splitting, Scenario(I)) or a weak inflection (weak splitting, Scenario(II)) in the conductance.

Figure 3(b) explores the effect of a constant dc offset bias on both the shape and position of the conductance feature. At $`V_{SD}`$ = 0 (black curves) we observe evolution of the feature from $`0.75`$ towards $`0.5\times 2e^2/h`$ with increasing $`V_T`$, consistent with the low temperature data for the $`l=0.5\mu m`$ quantum wire shown in Fig.(3a) The application of a DC offset bias $`V_{SD}`$ = 0.5mV causes the feature to rise towards $`0.75\times 2e^2/h`$, with a much weaker dependence on $`V_T`$. This behavior mirrors that in Fig 3(a) since increasing $`V_{SD}`$ or increasing $`T`$ will distribute electrons between both spin bands.

With the application of a dc source - drain bias the spin gap can be studied as a function 1D density, controlled by the side gate bias $`V_S`$. Figure 4 shows the differential conductance ($`di/dv`$) of a $`l=1\mu m`$ quantum wire as a function of $`V_{SD}`$ for different side gate voltages at $`T`$ = 50mK. We compare the differential conductance at two different 2D densities ($`V_T`$ = 385mV, $`n_{2D}=2.6\times 10^{11}cm^2`$ in (a) and $`V_T`$ = 700mV, $`n_{2D}=5.4\times 10^{11}cm^2`$ in (b)). In these plots conductance plateaus appear as a grouping of individual curves, as can clearly be seen in Fig. 4(b) where plateaus occur at $`2`$, $`1`$ and $`0.75\times 2e^2/h`$ for $`V_{SD}`$ = 0. In both Figs. 4(a) and (b) clear half-plateaus at $`0.5`$ & $`1.5\times 2e^2/h`$ can be seen developing near $`V_{SD}=\pm 1.5mV`$, (although the $`0.5`$ half-plateau on the right side of each graph is suppressed due to the asymmetric bias across the constriction near pinch-off moreno ).

Most importantly there are no plateaus near $`0.25`$ and $`1.25\times 2e^2/h`$, despite the presence of strong features at $`0.750.9\times 2e^2/h`$. This cannot be explained by a static spin gap, but is a natural consequence of a density dependent gap: at low 1D densities (small conductances) the spin gap has not yet developed, but at larger densities the gap opens up and features are observed at $`0.75\times 2e^2/h`$.

Further evidence for the density dependence of the gap can be seen in the region close to zero bias where a characteristic ‘cusp’ feature is seen below the $`1\times 2e^2/h`$ plateau. Traces associated with this feature start at $`0.5\times 2e^2/h`$ at zero source - drain bias (scenario I) and move towards $`0.75\times 2e^2/h`$ as $`V_{SD}`$ is increased. As the 1D density is increased, (by altering $`V_S`$) the spin-gap widens, and a larger $`V_{SD}`$ must be applied before the conductance increases above $`0.5\times 2e^2/h`$. The cusp feature is a result of many of these traces overlapping, and the strength and width of the feature is a measure of how large the spin gap is, and how rapidly it changes with 1D density. The increasing strength of the spin splitting with increasing $`V_T`$ can be seen in Figs. 4(a) and (b), where the cusp becomes deeper and more strongly resolved. With the application of an external parallel magnetic field the spin gap opens further and the cusp widens (Fig.4(c)) contacts . In future it may be possible to analyse the detailed shape of the cusp feature in order to obtain quantitative information about how the spin gap evolves with 1D density and parallel magnetic field.

In conclusion, we have presented a simple model to explain the $`0.7\times 2e^2/h`$ conductance feature in terms of a density dependent spin polarisation arising in the 1D region. While our phenomenological model is consistent both with experimental data for ultra low-disorder quantum wires presented here, and with other published data, a detailed microscopic explanation of the spin polarisation is still lacking. In particular, the rate at which the spin-splitting grows with 1D density is sample dependent and seems to depend on the length of the 1D region, the surface gate geometry, and the 2D electron density. A detailed understanding of the spin polarisation will have implications for 1D electron transport in mesoscopic devices and may have important applications to the field of spintronics.

We thank D. Barber, R. P. Starrett, and N. E. Lumpkin for technical support. This work was funded by the Australian Research Council.


text:
# Detection of the BL Lac Object H1426+428 at TeV Gamma Ray Energies

## 1 Introduction

Blazars are the main class of Active Galactic Nuclei (AGN) detected above 100 MeV by the EGRET experiment on the Compton Gamma Ray Observatory (CGRO) and by ground-based $`\gamma `$-ray observatories (Mukherjee et al., 1997; Weekes, 2001). They comprise a subclass of AGN and are characterized by a highly variable non-thermal continuum, strong variable optical polarization, the lack of a UV-excess (or “blue bump”) and a core-dominated radio morphology. BL Lacertae (BL Lac) objects are a subclass of blazars that are notable for their lack of prominent emission lines. The broad-band double-humped Spectral Energy Distributions (SEDs) of BL Lacs identified in x-ray surveys differ significantly from the SEDs of those identified in radio surveys. This led to the sub-classification of BL Lacs into High-frequency-peaked BL Lacs (HBLs) and Low-frequency-peaked BL Lacs (LBLs) based on the ratio of their x-ray to radio flux densities (Padovani & Giommi, 1995). The first “hump”, generally assumed to be the peak of the synchrotron emission (in a $`\nu F_\nu `$ representation), is in the IR-optical for LBLs and in the EUV-soft x-ray band for HBLs. BL Lac objects make up a significant fraction of the 70 blazars in the 3rd EGRET Catalog (Hartman et al., 1999) and most of them are classified as LBLs. It has been shown that there is not a sharp division between these two classes of objects (see, e.g., Fossati et al., 1998; Ghisellini, 1999).

Recent studies with ground-based $`\gamma `$-ray telescopes have produced evidence for TeV $`\gamma `$-ray emission from seven BL Lac objects, five of which are classified as HBL and two as LBL (Weekes, 2001). The most prominent of these are Markarian 421 (Mrk 421), which has been detected by five ground based imaging atmospheric Čerenkov $`\gamma `$-ray observatories, and Markarian 501 (Mrk 501), which has been detected by six such observatories. The emission from both of these objects can be explained by Compton-synchrotron models although detailed modeling is still fraught with many uncertainties. Both objects are characterized by rapid variability on time-scales from hours to months. In the TeV energy range, the energy spectrum of Mrk 501 and, more recently, that of Mrk 421, have been shown to exhibit absorption-like features (Krennrich et al., 1999; Aharonian et al., 1999; Krennrich et al., 2001). The temporal and spectral properties of the other TeV BL Lacs are less well defined.

Since 1992, the Whipple Gamma Ray collaboration, using the 10m imaging atmospheric Čerenkov telescope on Mt. Hopkins, has been searching for TeV $`\gamma `$-ray emission from AGN. Initially the search was concentrated on blazars detected by EGRET at any redshift; these observations led to the detection of Mrk 421 (Punch et al., 1992) and upper limits on some 30 other blazars (Kerrick et al., 1995). More recently, the search has concentrated on nearby BL Lacs leading to the detection of Mrk 501 (Quinn et al., 1996) and 1ES 2344+514 (Catanese et al., 1998). Between 1995 and 1998, the survey included 24 objects (17 HBLs, 7 LBLs) ranging in redshift from 0.046 to 0.44; the results for these objects will be published shortly (D. Horan et al., in preparation). Although none of these observations resulted in a detection, the observations of H1426+428 yielded the highest consistently positive statistical significances.

In this paper, evidence is presented for the detection of very high energy (VHE) $`\gamma `$-rays from H1426+428. The H1426+428 observations are divided into two categories - those taken as part of the general blazar survey between 1995 and 1998, and the subsequent concentrated observations carried out between 1999 and 2001. A preliminary energy spectrum is derived and the implications of the detection of H1426+428 at these energies are discussed.

## 2 H1426+428 at Other Wavelengths

H1426+428 was discovered in the 2-6 keV band by HEAO-1 (Wood et al., 1984) and was classified as a BL Lac object in 1989 (Remillard et al., 1989). It has a cosmological redshift of 0.129 and is located in the constellation of Böotes ($`\alpha _{J2000}`$ = 14h 28m 32s.7, $`\delta _{J2000}`$ = +42 40 20′′). It is an optically faint object (mV = 16.9) and is believed to be at the center of an elliptical galaxy. H1426+428 is bright in the x-ray band with a 2-6 keV luminosity of $``$ 1044 erg s-1 which is typical of the BL Lacs found with HEAO-1. Both the flux of H1426+428 in the 2-10 keV band and its spectral index above 2 keV have been found to change over time (Costamante et al., 2001).

Three recent BL Lac surveys, DXRBS (Perlman et al., 1998), RGB (Laurent-Muehleisen et al., 1999), and REX (Caccianiga et al., 1999), have shown that, despite the conventional subdivision into HBLs and LBLs, BL Lacs actually form a continuous class with respect to the peak of the synchrotron emission, which smoothly ranges between IR and soft x-ray frequencies and up to the 2-10 keV band for sources like Mrk 421. Beppo-SAX observations of Mrk 501 (Pian et al., 1998) and 1ES 2344+514 (Giommi, Padovani, & Perlman, 2000) revealed that, at least in a flaring state, the first peak of the SED can reach even higher energies, at or above 100 keV.

In 1998-1999, Beppo-SAX performed an observing campaign with the aim of finding and studying other sources as “extreme” as Mrk 501 is in its flaring state. The candidates for the Beppo-SAX survey were selected from the Einstein Slew Survey and the RASSBSC catalogs. These Beppo-SAX observations (Costamante et al. 2000a, 2000b, 2001) revealed four new “extreme” HBLs, selected because they have high synchrotron peak frequencies and are therefore possible TeV emitters. These four candidates for TeV emission were: 1ES 0120+340, PKS 0548-322, 1ES 1426+42.8 (i.e. H1426+428) and H2356-309.

The spectra for three of these objects (1ES 0120+340, PKS 0548-322, and H2356-309) were well fitted by a convex broken power law with a break energy, and hence the peak of the synchrotron emission, occurring at about 1.4 keV for 1ES0120+340, 4.4 keV for PKS0548-322 and 1.8 keV for H2356-309. In the case of H1426+428 however, no evidence for a spectral break up to 100 keV was found. Instead, its spectrum was well fitted by a single power-law, with a flat spectral index of 0.92 up to 100 keV, thus constraining the peak of the synchrotron emission to lie near or above this value. The best fit of a pure homogeneous synchrotron self-Compton (SSC) model for H1426+428 (Costamante et al., 2000a) predicted detectable $`\gamma `$-ray emission at TeV energies. At the time of the Beppo-SAX observation the observed x-ray flux in the 2-10 keV band was at one of the lowest levels ever recorded from H1426+428, indicating that it was not in a flaring state (Costamante et al., 2001). This implies that, in the event of a flare, the synchrotron peak could shift to even higher values as was observed for both Mrk 501 and 1ES 2344+514. It also indicates that highly relativistic electrons are present which makes H1426+428 the most promising candidate for TeV emission from this survey.

## 3 Observing Technique

The observations reported in this paper were taken with the 10m reflector at the Whipple Observatory on Mount Hopkins in southern Arizona (elevation 2.3 km) using the atmospheric Čerenkov imaging technique. During the course of the observations presented here, many changes were implemented on the 10m telescope. The imaging camera (Cawley et al., 1990) was upgraded a number of times (Finley et al., 2000), the number of photo-multiplier tubes (PMTs) being increased with each iteration, resulting in the current high resolution, 490 pixel camera (Finley et al., 2001). The triggering electronics were upgraded to a Pattern Selection Trigger (Bradbury et al., 1999) and light cones were installed in front of the PMTs to increase their light collection efficiency. The unprotected, anodized, front-aluminized mirrors on the 10m were recoated during the timespan of the observations. Because of changes in PMT configuration, triggering, mirror reflectivity, and light collection efficiency, the energy response of the camera for $`\gamma `$-ray detection varied. Since various changes were made at the beginning of each observing season (typically between June and September), it is convenient to consider each observing season separately and to characterize it by the observed response from the assumed standard candle, the Crab Nebula. The different configurations of the camera and the resultant peak response energies are summarized in Table 1.

The Čerenkov light images from each shower are analyzed off-line using a moment analysis (Reynolds et al., 1993). The derived image parameters are used to distinguish candidate $`\gamma `$-rays from the large background of cosmic-rays. These parameters include length, width, distance (from the optic axis), and orientation (alpha). In addition, the two highest recorded signals in individual pixels (sig1 and sig2) are noted as well as the total light in the image (size). Monte Carlo simulations were used to determine the approximate limits on the parameters to be used for the identification of candidate $`\gamma `$-ray images; these limits were then optimized on contemporaneous observations of the Crab Nebula. The results of an analysis are graphically presented as a histogram of the alpha parameter - usually referred to as an alpha plot. In such histograms, the alpha values for all the events that passed all $`\gamma `$-ray selection criteria except for the alpha cut, are plotted. For a $`\gamma `$-ray source, an excess would be expected at low values of alpha.

Two modes of observation were used: on/off and tracking. In on/off mode, a 28 sidereal minute run is taken with the candidate $`\gamma `$-ray source at the center of the field of view - the on run. The off run, a ’control’ run, is also taken for 28 sidereal minutes through the same range of azimuth and elevation, thus enabling this region of the atmosphere to be characterised in the absence of the $`\gamma `$-ray candidate. In the tracking mode, the object is continuously tracked for 28 sidereal minutes with no control run being taken.

Even in the absence of a $`\gamma `$-ray source, a certain percentage of events recorded will pass all of the $`\gamma `$-ray selection criteria. This background level of events, which depends on a number of factors including sky brightness, elevation, and weather, needs to be established in order to calculate the statistical significance of any apparent excess of $`\gamma `$-rays detected. There are different methods to estimate this background depending on which mode of observation was used.

For data taken in on/off mode, the off scan provides an estimate of the number of $`\gamma `$-ray like events that would be recorded from this range of azimuth and elevation, in the absence of the candidate $`\gamma `$-ray source. In this data acquisition mode, differences in night-sky brightness between the on and off regions of the sky can introduce a bias when the data are analyzed. In order to reduce this bias, the standard deviations of the night-sky background in each PMT are compared for the on and the off runs. Gaussian noise is then added in quadrature to the signal from whichever region of the sky (on or off) is the darker, so as to match the standard deviations for the tube in the on and off runs. This technique is known as software padding and is described in Cawley (1993).

Unlike observations taken in the on/off mode, scans taken in the tracking mode do not have independent control data. A background estimate is essential in order to predict what number of the background events, passing all cuts, would have been detected during the scan in the absence of the candidate $`\gamma `$-ray source. Since most of the H1426+428 observations were taken in the tracking mode, two different methods of background estimation were used. Each of these methods is outlined below.

In the first method, which is the standard analysis method used for data taken in the tracking mode (Catanese et al., 1998), the alpha plot was characterized when there was no $`\gamma `$-ray source at the center of the field of view. This was done by analyzing ’darkfield data’ which consisted of data taken in the off mode, and of observations of objects found not to be sources of $`\gamma `$-rays. A large database of these scans were analyzed, giving the shape that an alpha plot would be expected to have when no $`\gamma `$-rays are present.

Most $`\gamma `$-rays from an object at the center of the field of view will have small values of the *alpha* parameter. Hence, the *alpha* distribution beyond 20 can be assumed to be independent of the $`\gamma `$-ray source, and thus representative of the background level of $`\gamma `$-ray-like events in the field of view. Using the darkfield data, a ratio was calculated to scale the number of events between 20 and 65 to the number that pass the *alpha* cut. This ratio, the ’tracking ratio’ ($`\rho `$), was then used to scale the 20 to 65 region of the alpha plot for the tracking scan, to estimate the background level of events passing all cuts. Its value and statistical uncertainty ($`\mathrm{\Delta }\rho `$) are given by:

$$\rho \pm \mathrm{\Delta }\rho =\frac{N_{alpha}}{N_{control}}\pm \sqrt{\frac{N_{alpha}}{N_{control}^2}+\frac{N_{alpha}^2}{N_{control}^3}}$$
(1)
where $`N_{alpha}`$ is defined as the number of events in the darkfield data that pass all the $`\gamma `$-ray selection criteria including the *alpha* cut, while $`N_{control}`$ is the number of such events with *alpha* between 20 and 65. The tracking ratio was calculated using the darkfield data available for each observing year over the same range of elevation angles as the source observations. In addition, the tracking ratio was checked using the off data on H1426+428 from on/off runs (where available); within statistics these tracking ratios were consistent with the standard tracking ratio derived from the full yearly database. To check for systematic effects associated with the determination of the global tracking ratio, the darkfield data were subdivided on the basis of the region of the sky which they were from. Each of these dark-sky regions was checked individually and the tracking ratio determined and compared to the global tracking ratio. It was found that the global tracking ratio and these individual tracking ratios were consistent within statistical uncertainty.

In the alternative method for estimating the background, each H1426+428 tracking scan was assigned a suitable off scan as its background. The off database comprised all off scans taken during the observing season including H1426+428 off runs. In order to establish the most suitable off scan for each H1426+428 tracking scan, each off run was characterized by parameters pertinent to the conditions during the scan. These consisted of the elevation, the number of PMTs switched off and the mean night-sky background during that scan, along with the throughput factor recorded on the night of the scan and the date on which the scan was taken. The throughput factor (Holder et al., 2001) is a measure of sky clarity based on the spectrum of the total amount of light produced by background cosmic rays. The same information was assembled for each H1426+428 tracking scan. Then, for each H1426+428 tracking scan, the off scan whose conditions most closely matched it was deemed the most suitable background estimate, and was used as the background for that H1426+428 tracking run. In cases where there was not a suitable off run, the H1426+428 tracking run was omitted. With each H1426+428 tracking scan then having an associated off run, the data were analyzed as if they were taken in the on/off mode, and hence, software padding was applied. The sky quality changes from night to night, causing variations of $`10\%`$ in the raw trigger rate, for data taken on the same source at the same elevation. Therefore, even after selecting off scans with run conditions which match those of the on scans, it is still necessary to normalize the total number of events on and off . We calculate the normalizing factor from the ratio of the number of events in the 30 \- 90 control region of the on and off alpha plots.

It was found that both methods of background estimation were consistent, and indicated excesses of similar significance, thus suggesting that the $`\gamma `$-ray excess can be reliably determined for tracking scans.

## 4 Database

For the purpose of data analysis, two different strategies were employed: (a) All of the data taken in the on and the tracking modes were combined and subjected to a tracking analysis. (b) All of the data taken in the on mode, along with their corresponding off data, were analyzed together using the pairs analysis. The $`\gamma `$-ray selection criteria, including the alpha bound, were re-optimized each year using data from the Crab Nebula. The optimum value for the alpha bound was found empirically to increase from 10 to 15 when the smaller, high resolution 490 pixel camera was installed in the summer of 1999.

### 4.1 1995 - 1998 Observations

In June-July, 1995, H1426+428 was observed for 3.5 hours in the tracking mode with the 109 pixel camera (spacing 0.259, 3.0 field of view and standard trigger). In February-June, 1997, another 13.2 hours were used to observe H1426+428 in the tracking mode; the camera configuration was as before but the number of pixels was increased to 151. In April, 1998, 0.87 hours of H1426+428 data were taken; the 331 pixel camera was installed at this time with the standard trigger.

### 4.2 1999 Observations

During the 1999 observing season, a total of 7 on/off pairs and 51 tracking runs were taken with the 2-fold Pattern Selection Trigger. Three tracking runs were excluded from further analysis based on fluctuations in the raw rate and inferior weather conditions. This left a total of 24.35 hours (7 on runs and 48 tracking runs) of data available for analysis. These observations were initiated as part of the general BL Lac survey (D. Horan et al., in preparation) with extra data being taken in 1999 because of an initial indication of a signal from H1426+428. A tracking ratio of 0.232 $`\pm `$ 0.005 was determined for *alpha* $`<`$ 10.

### 4.3 2000 Observations

During the 2000 observing season, 33 on/off pairs and 35 tracking runs were taken on H1426+428. Eleven on source runs were excluded from further analysis based on fluctuations in the raw rate and inferior weather conditions, leaving a total of 26.46 hours of data available for tracking analysis (57 tracking/on) and 13.86 hours of data in the on/off mode (30 on/off pairs). A tracking ratio of 0.312 $`\pm `$ 0.002 for alpha $`<`$ 15 was determined. H1426+428 off data were analyzed to ensure that this tracking ratio was consistent with the tracking ratio derived for these data.

### 4.4 2001 Observations

A total of 42.24 hours were spent on the source during 2001; this comprised 39 on/off pairs and 59 tracking runs. After excluding data taken at low elevation and in unsatisfactory weather conditions, 38.10 hours of data remained (34 on/off pairs and 53 tracking runs). The tracking ratio for this observing season, again for alpha $`<`$ 15, was calculated to be 0.323 $`\pm `$ 0.002. As with the 2000 data, checks were performed to ensure that this ratio was appropriate for the H1426+428 data. A summary of the data taken and analyzed in the 2000 and 2001 observing seasons is give in Table 2.

## 5 Results

### 5.1 1995 - 1998 Observations

The observations of H1426+428 between 1995 and 1998 were analyzed as part of a general BL Lac survey. The results are summarized in Table 3. Although not statistically significant, H1426+428 was one of only two objects, out of the 24 objects observed in this survey, to show consistently positive results.

### 5.2 1999 Observations

A total of 24.35 hours of data were taken on H1426+428 during 1999. The net excess of +0.9$`\sigma `$ from the combined 1999 observations was not statistically significant.

### 5.3 2000 Observations

A tracking analysis of all the on-source data taken on H1426+428 during 2000, using the standard tracking ratio calculated for that observing season, gave a +4.2$`\sigma `$ excess corresponding to a $`\gamma `$-ray rate of 0.24 $`\pm `$ 0.06 $`\gamma `$ min-1. When these data were analyzed independently using the alternative method of background estimation, the excess was at the +3.1$`\sigma `$ level with a $`\gamma `$-ray rate of 0.21 $`\pm `$ 0.07 $`\gamma `$ min-1. The alpha plot for these data is shown in Figure 1 (left panel) along with the alpha plot from the matched off data. The background shown on the figure (dashed line) consists of the off data which were chosen to match the characteristics of the H1426+428 tracking data. These data are scaled to match the H1426+428 data in the 30 \- 90 region of the alpha plot. The net significance for the 30 on/off pairs taken in 2000 was +1.2$`\sigma `$, with a $`\gamma `$-ray rate of 0.10 $`\pm `$ 0.09 $`\gamma `$ min-1. The alpha plot for these data is shown in Figure 1 (right panel). On the night of May 30, 2000 (MJD 51694), the first tracking scan on H1426+428 gave a significance of +3.0$`\sigma `$, using the standard tracking ratio method of background estimation. This prompted observers to take four more tracking scans which led to the detection of a signal at the +3.7$`\sigma `$ for this night. The $`\gamma `$-ray rate for these 5 scans (2.15 hours) was 0.67 $`\pm `$ 0.18 $`\gamma `$ min-1.

### 5.4 2001 Observations

A tracking analysis of all of the on-source data taken on H1426+428 during 2001, using the tracking ratio calculated for this observing season, gave an excess of +5.4$`\sigma `$ corresponding to a $`\gamma `$-ray rate of 0.36 $`\pm `$ 0.07 $`\gamma `$ min-1. An independent analysis of 31.13 hours of these data, using the alternative method of background estimation, resulted in a +5.5$`\sigma `$ excess with a $`\gamma `$-ray rate of 0.44 $`\pm `$ 0.08 $`\gamma `$ min-1. The alpha plot for these data is shown in Figure 2 (left panel) along with the alpha plot from the matched off data. As before, these data are scaled to the H1426+428 data between alpha values of 30 and 90.

The significance from the 34 on/off pairs taken during the 2001 observing season is +2.0$`\sigma `$, with a $`\gamma `$-ray rate of 0.22 $`\pm `$ 0.11 $`\gamma `$ min-1. The alpha plot for these data is shown in Figure 2 (right panel). A summary of the results of the observations made on H1426+428 during 2000 and 2001, using the alternative method of background estimation, is presented in Table 4.

When both the differences in exposure time and methods of background estimation are taken into account, the results from the on/off pairs and the tracking data taken on H1426+428 during each observing season are consistent, at the 2$`\sigma `$ level, with each other.

### 5.5 Comparison of the Gamma-Ray and X-Ray Flux from H1426+428

The $`\gamma `$-ray rates for the 2000 and 2001 H1426+428 data were compared with the x-ray flux from the All Sky Monitor (ASM) instrument on board the Rossi X-ray Timing Explorer (RXTE; Levine et al., 1996). A correlation was sought between the nightly $`\gamma `$-ray rates and the “one-day average” x-ray data points. Only nights on which there were both x-ray and $`\gamma `$-ray data were considered when performing the correlation. This left 37 nights for analysis during 2000, and 39 nights during 2001.

Since $`\gamma `$-ray observations can only be taken on moonless nights, there are a few ($``$ 6) nights around the time of full moon each month, when no observations can be made. The periods during which the $`\gamma `$-ray data are taken are referred to as ’darkruns’. Correlations between the x-ray and $`\gamma `$-ray rates were sought both over the entire observing season, and during darkruns with more than two nights on which both x-ray and $`\gamma `$-ray data were taken. No evidence for significant nightly correlation was found for either observing season. The x-ray and $`\gamma `$-ray rate curves are shown in Figure 3. The rates plotted here are the average monthly rates for the 2000 and 2001 data. The $`\gamma `$-ray rates shown were calculated using the alternative method of background estimation. There is some evidence for an anti-correlation between these average monthly rates, especially in the 2001 data.

## 6 Spectral Characteristics of H1426+428

The TeV flux from H1426+428 is weak and near the threshold of sensitivity of the 10m $`\gamma `$-ray telescope. The photon flux is so small that it is impractical to apply standard spectral analysis techniques, for example, Mohanty et al. (1998), to these data. In fact, only the differential flux from the source for the energy at which the telescope is most efficient at detecting $`\gamma `$-rays, can be evaluated with reasonable accuracy. Such a flux estimate accounts mostly for the photons with energies corresponding to the peak in the differential detection rate of the telescope. This energy, Ep, depends on the spectral index of the source in question. For observations of H1426+428 with the Whipple telescope in 2001, $`E_p`$ was found to be between 280 - 360 GeV. The rate of change of the number of excess events, presumably photons from the source, with the total amount of light in the Čerenkov image, the size, is directly related to the spectral index of the source. This relationship is established using Monte Carlo simulations; the spectral index of the source in the vicinity of $`E_p`$ is then evaluated. This in turn allows estimates of $`E_p`$ for this spectral index, and the differential and integral fluxes to be improved by reducing the uncertainty due to the unknown spectral index.

Figure 4 shows the excess events detected from the direction of H1426+428 as a function of integrated Čerenkov light in the shower image. For comparison, this plot also shows the excess events detected from the Crab Nebula after a 4.1 hour exposure, which was chosen to produce a similar maximal excess of $``$ 103 events as was found during the 32.5 hours of observation of H1426+428 during 2001. The lowest size cut of 366 digital counts is the result of signal to noise optimization of the H1426+428 data, producing an overall significance close to +6$`\sigma `$. The highest size cut, 840 digital counts, is limited by the rapidly declining significance to the level of +3.5$`\sigma `$. To a first approximation, the dependence of the excess events on the logarithm of the size cut applied, is a linear function whose slope is monotonically related to the spectral index of the source. This plot indicates that the H1426+428 spectrum is softer than the well-studied spectrum of the Crab Nebula.

Monte Carlo simulations were used to relate the observed parameters of linear fits to the data in Figure 4, to the spectral characteristics of H1426+428, and the Crab Nebula. The statistical correlation of the errors shown in this integral plot of excess events has been taken into account to derive optimal estimates of the spectral indices, differential and integral fluxes, and their errors. A summary of the results is given in Table 5. For the Crab Nebula a spectral index of 2.75 with an uncertainty of $`6`$% (1$`\sigma `$) was derived. The peak of the differential detection rate of the Whipple telescope from a source with such a spectral index has been found to be around 360 GeV (2001 observing season). Based on the small Crab Nebula data-set, the spectral index as well as the differential and integral fluxes have been derived, and are found to be consistent, at the 2$`\sigma `$ level, with our previous observations (Hillas et al., 1998). Analysis of the H1426+428 data-set indicates a substantially steeper spectral index of 3.55 with a 13% relative error. For this spectral index and the 366 digital counts size cut applied in the data analysis, the differential photon detection rate peaks at $``$ 280 GeV. Table 5 shows estimates of the fluxes at 390 GeV assuming the same spectral index. This energy was chosen by increasing the size cut until the signal to noise ratio decreased to $``$ 4.5 $`\sigma `$; below this value the signal was not considered strong enough for reliable calculation of the H1426+428 flux.

The observed energy fluxes of H1426+428 at $`280`$ and $`390`$ GeV are shown in Figure 5 together with the predicted SED for this source suggested by Costamante et al. (2001). This figure also shows the tentative detection (+3.1 $`\sigma `$) made by the Whipple 10m telescope during 2000. Operating at a substantially higher peak response energy at that time, the Whipple telescope observed this source at Epeak $``$ 430 GeV. Based on these data an estimate of the energy flux of 0.86 $`\pm `$ 0.33 $`\times `$ 10-11 erg cm-2 s-1 has been derived. Due to the very weak signal however, we have not attempted to find the spectral index of H1426+428 for this data-set. As a result, the indicated error is only statistical and it does not include errors associated with the unknown spectral properties of the source. The $`430`$ GeV peak energy itself is subject to a similar uncertainty.

## 7 Discussion

Observations of H1426+428 with the Whipple 10m telescope since 1995 have returned consistently positive excesses although the statistical significance has not always been high. H1426+428 was on the very short list of predicted TeV emitters published by the *Beppo*SAX collaboration (Costamante et al., 2000a) and was also singled out as the most probable TeV emitter amongst these due to the high frequency of its synchrotron peak. It was observed as part of a BL Lac survey at the Whipple Observatory between 1995 and 1998, and was singled out for more extensive observations between 1999 and 2001. The strongest evidence for VHE emission was recorded during the 2000 and 2001 observing seasons, with signals at the 3.1$`\sigma `$ and 5.5$`\sigma `$ levels, respectively, being detected. Although it is difficult to combine results from different configurations of the camera, the combined significance of the results from 2000 and 2001 is clearly equivalent to a greater than 5$`\sigma `$ detection, the usual level required to establish the existence of a new TeV source.

There have been two occasions (in 1999 and in 2000) when there was marginal evidence for a transient signal from H1426+428. In 1999, after 2.3 hours of observations (over two nights) a signal at the +3.1$`\sigma `$ level was recorded while in 2000, a signal at the +3.7$`\sigma `$ level was recorded after 2.2 hours of observations (on a single night).

After 18.1 hours of observations on H1426+428, the CAT collaboration reported a 3$`\sigma `$ upper limit of 3.54 x 10-11 cm-2 s-1 at an energy threshold of 250 GeV (Piron, 2000). These observations were carried out at the beginning of February 1999 in response to an announcement from the *Beppo*SAX and the RXTE/ASM teams, which revealed H1426+428 to be transitioning to a high state.

If independently verified, the observations reported here represent an important addition to the catalog of TeV-emitting BL Lacs (Table 6). That the source was predicted to be a TeV $`\gamma `$-ray source based on its x-ray spectrum is important in that it signifies the maturity of the observational techniques and the theoretical understanding of BL Lacs. It reinforces the symbiosis between observations at x-ray wavelengths, particularly in hard x-rays, and those at TeV energies, particularly those with good sensitivity below energies of 1 TeV. The existence of a population of sources whose most prominent emission is at energies of 10-100 keV and 300-1000 GeV points to a fruitful overlap between the next generation of ground-based atmospheric Čerenkov telescopes (Weekes et al., 2001) and the future hard x-ray experiment, EXIST (Grindlay et al., 2000). It is remarkable also in that this source, like Mrk 501 and 1E2344+514, is not included in the 3rd EGRET Catalog, supporting the interpretation of Ghisellini (1999) that in blazars sufficiently powerful to be detected readily by EGRET, the spectrum will terminate at lower energies and not extend to the TeV domain.

Although the properties of H1426+428 reported here are scanty in comparison with the other, better studied TeV BL Lacs, they agree in principle with the characteristics now well established for Mrk 421 and Mrk 501 (and to a lesser extent, 1ES 2344+514). Although the flux from H1426+428 is low, there is some evidence for time variability. There is also an indication of a hard $`\gamma `$-ray spectrum. There is no evidence for a correlation of x-ray and TeV $`\gamma `$-ray fluxes (unlike Mrk 501) but the sensitivity of observations to date is limited and even in the stronger sources this correlation has been shown to be complicated.

Ghisellini (1999) has proposed that there is a continuous sequence of BL Lacs with the peak of the synchrotron spectrum increasing as the luminosity decreases. H1426+428 is another example of an ”extreme” BL Lac characterized by a synchrotron spectrum that may peak near 100 keV and a relatively weak luminosity. These objects are the best candidates for TeV emission based on Compton-synchrotron models of BL Lac jets. However, the observed variations in sources such as Mrk 501 cannot be accounted for by simple one-zone homogeneous self-Compton models and may require another source of optical target photons. It will require the detection of a greater population of sources such as H1426+428 to test these models fully. It has been noted (Ghisellini, 1998) that these TeV-emitting BL Lacs are also characterized by a relatively strong radio luminosity and that the radio luminosity may be a measure of the density of additional seed photons at infrared-optical wavelengths needed to explain the TeV emission.

This is the most distant of the TeV-detected BL Lacs classified as HBL, and hence it has promising implications for the detection of more BL Lacs at $`z`$ $`>`$ 0.1. Stronger detections of such sources, which allow an accurate measure of the TeV energy spectrum, may place significant limits on the density of the intergalactic background light. The steep spectrum derived here is consistent with many models of intergalactic absorption but could also be intrinsic to the source.

Interpretation of the theoretical prediction of the intrinsic spectrum of H1426+428 and the observational results reported in this paper are complicated by the possible strong attenuation of the high energy photons by the diffuse intergalactic infrared background. The appearance of such a cutoff in the spectra of TeV extragalactic sources was suggested in one of the pioneering works on TeV $`\gamma `$-ray absorption, by Stecker, de Jager & Salamon (1992), for the quasar 3C 279 ($`z`$ = 0.54). A later paper (V. V. Vassiliev et al., in preparation) will discuss this subject in detail in relation to H1426+428; here we note only that at 280 GeV the attenuation optical depth can be anywhere between $`z`$ of 0.1 and 1.0 depending on the extreme upper and lower limits known for the density of the extragalactic background light. Due to this potentially large absorption effect it is possible, for example, that the intrinsic spectral index of this source is softened by $``$ 1.75 to produce the observed value of 3.55, if the intrinsic properties of H1426+428 are similar to those of Mrk 501 (Vassiliev et al., 2001). It remains certain however, that being the most distant HBL detected at sub-TeV energies, the observed spectral properties of H1426+428 and the understanding of its intrinsic characteristics from multiwavelength observations will provide the most constraining data for extragalactic background light studies.

## 8 Acknowledgments

The authors would like to thank Kevin Harris, Joe Melnick, Emmet Roach, and all the staff at the Whipple Observatory for their support. We also thank the referee for many constructive and useful suggestions. This research was supported in part by the U. S. Department of Energy, PPARC, and Enterprise Ireland.

## Appendix A Data Log for 2000 and 2001 Observations