text:
# Lyapunov exponents in continuum Bernoulli-Anderson models

## 1. The Main Result

We study Anderson-type random Schrödinger operators

$$H_\omega =\frac{d^2}{dx^2}+\underset{n}{}q_n(\omega )f(xn)$$
(1)
in $`L^2()`$. The single site potential $`f`$ is assumed to be real, supported in $`[1/2,1/2]`$, locally in $`L^1`$, and not identical to $`0`$ (in $`L^1`$-sense). The coupling constants $`q_n(\omega )`$, $`n`$, are independent, identically distributed Bernoulli random variables, i.e. they have distribution $`\mu `$ with supp$`\mu =\{0,1\}`$. Thus the random potential $`V_\omega (x)=_nq_n(\omega )f(xn)`$ is in $`L_{loc,unif}^1()`$ for all $`\omega `$, which allows to associate a unique selfadjoint operator with $`H_\omega `$ that may equivalently be defined by form methods or in the sense of Sturm-Liouville theory.

Our interest in Bernoulli-Anderson models of this type arises from the fact that their spectral properties are not as well understood, by mathematically rigorous standards, as those of Anderson models with continuous (or absolutely continuous) distribution $`\mu `$. Furthermore, these discrete distributions are physically more relevant as they are modeling the charge numbers of nuclei.

One of the central objects in the study of one dimensional random operators is the Lyapunov exponent. To define it for our model, let $`g_E(n,\omega )`$ be the transfer matrix of

$$u^{\prime \prime }+V_\omega u=Eu$$
(2)
from $`n1/2`$ to $`n+1/2`$, i.e. for any solution of (2) one has

$$\left(\begin{array}{c}u(n+1/2)\\ u^{}(n+1/2)\end{array}\right)=g_E(n,\omega )\left(\begin{array}{c}u(n1/2)\\ u^{}(n1/2)\end{array}\right).$$

Also, for $`n`$ let $`U_E(n,\omega )=g_E(n,\omega )\mathrm{}g_E(1,\omega )`$, and define the Lyapunov exponent at $`E`$ by

$$\gamma (E)=\underset{n\mathrm{}}{lim}\frac{1}{n}𝔼(\mathrm{log}U_E(n,\omega )),$$
(3)
where $`𝔼`$ denotes expectation with respect to $`\omega `$. Existence of $`\gamma (E)`$ follows from the subadditive ergodic theorem, e.g. . Since $`U_E1`$, we have $`\gamma (E)0`$, and it can also be seen that one gets the same $`\gamma (E)`$ if transfer matrices on the negative half line are used analogously in (3).

For many applications it is crucial to know at which energies one has that $`\gamma (E)>0`$, which corresponds to exponential growth (or decay) of the solutions to (2). Our main goal here is to show that under the above assumptions this holds for all but a discrete set $`M`$ of energies $`E`$, and to explicitly describe $`M`$ in terms of the transmission and reflection coefficients for scattering at the single site potential $`f`$. To define them, let $`k\{0\}`$ and $`u_+(x,k)`$ be the Jost solution of

$$u^{\prime \prime }+fu=k^2u,$$
(4)
i.e. the solution satisfying

$$u_+(x,k)=\{\begin{array}{cc}e^{ikx}\hfill & \text{for }x1/2\text{,}\hfill \\ a(k)e^{ikx}+b(k)e^{ikx}\hfill & \text{for }x1/2\text{.}\hfill \end{array}$$
(5)
Since $`e^{ikx}`$ and $`e^{ikx}`$ are linearly independent solutions of $`u^{\prime \prime }=k^2u`$, this defines $`a(k)`$ and $`b(k)`$ uniquely. They are related to the transmission and reflection coefficients used in physics by $`a=t^1`$ and $`b=rt^1`$. In particular, vanishing of $`b`$ is equivalent to vanishing of $`r`$. For real $`k`$ we use constancy of the Wronskian to get

$$|a(k)|^2|b(k)|^2=1,$$
(6)
corresponding to the familiar $`|r|^2+|t|^2=1`$.

It can also be seen from (5) that

$$\left(\begin{array}{c}a(k)\\ b(k)\end{array}\right)=\frac{1}{2ik}e^{ik/2}\left(\begin{array}{cc}ik& 1\\ ike^{ik}& e^{ik}\end{array}\right)\left(\begin{array}{c}u_+(1/2,k)\\ u_+^{}(1/2,k)\end{array}\right).$$

Since $`u_+`$ is a solution of (4), a linear differential equation which is analytic in $`k`$, and satisfies the analytic initial condition $`(u_+(1/2),u_+^{}(1/2))=(1,ik)`$, we see that $`u_+(1/2,k)`$ and $`u_+^{}(1/2,k)`$ are entire in $`k`$. Thus $`a(k)`$ and $`b(k)`$ are analytic in $`\{0\}`$ with a possible pole at $`k=0`$. Neither $`a(k)`$ nor $`b(k)`$ vanish identically. For $`a(k)`$ this is a trivial consequence of (6). For $`b(k)`$ it follows from the fact that a compactly supported $`f0`$ can not be reflectionless (i.e. a soliton), as follows from inverse scattering theory, e.g. , see also . Thus the roots of $`a(k)`$ and $`b(k)`$ can not accumulate, neither at (the pole) $`0`$ nor away from $`0`$.

We are now ready to state our main result:

###### Theorem 1.1.

Let $`EM`$, where the set of critical energies $`M`$ is given by

$`M`$ $`:=`$ $`\{\left({\displaystyle \frac{n\pi }{2}}\right)^2:n\}\{E=k^2:k>0,b(k)=0\}`$ (8)
$`\{\alpha ^2<0:a(i\alpha )a(i\alpha )b(i\alpha )b(i\alpha )=0\},`$
then $`\gamma (E)>0`$.

Note that the critical set $`M`$ is the union of the discrete set $`\{(n\pi /2)^2:n\}`$ and all real numbers in $`\{k^2:k\{0\},a(k)=0\text{or}b(k)=0\}`$ (this uses that for real $`k`$, $`a(k)0`$ and $`b(k)=\overline{b(k)}`$). Therefore $`M`$ is discrete.

Our main tool for proving Theorem 1.1 in Section 2 will be Fürstenberg’s Theorem, which has been used extensively in proofs of positivity of $`\gamma `$ for discrete one dimensional random operators, see \[12, Sec.14A\] for a summary. For continuum models, the use of Fürstenberg’s Theorem until recently seems to have been restricted to the special cases $`f=\chi _{[1/2,1/2]}`$ and $`f=\delta _0`$, a $`\delta `$-point-interaction . The first to have used Fürstenberg’s Theorem systematically for continuum Anderson models with general classes of single sites $`f`$ has been Kostrykin and Schrader . While they state their general results for absolutely continuous distribution $`\mu `$, they point out that most of their ideas can also be applied to discrete distributions. Our work here can be seen as an implementation of this fact with, as we feel, a minimal amount of technical effort.

A nice feature of using Fürstenberg’s Theorem in the proof of Theorem 1.1 is that once we know that the Theorem holds under the assumption that supp$`\mu =\{0,1\}`$, then it is immediately clear that the Theorem extends to the case where $`\{0,1\}`$ supp$`\mu `$. In fact, for larger support the critical set $`M`$ should be smaller. Kostrykin and Schrader note in  that $`M`$ should be empty if supp$`\mu `$ has at least one non-isolated point, but there doesn’t seem to be a proof of this yet. While Theorem 1.1 does not state that $`\gamma (E)=0`$ for all $`EM`$ (which in fact is not generally true), we will demonstrate in Section 3 that for Bernoulli-Anderson models one can indeed find many critical energies with $`\gamma (E)=0`$. In the discrete case this has already been observed for the so-called dimer model, see . In particular, Section 3 discusses the example $`f=\lambda \chi _{[1/2,1/2]}`$, $`\lambda `$ a constant, which leads to two different types of critical energies. They can be classified by the asymptotics of solutions of (2). In some cases they are of plane wave type, in particular bounded, while other critical energies lead to solutions which grow like $`\mathrm{exp}\sqrt{x}`$, due to a connection with random walks.

In  we use a result much like Theorem 1.1 to prove exponential and dynamical localization for continuum Bernoulli-Anderson models, where in addition we build on ideas which were developed for discrete models in . Discreteness of the critical set $`M`$ is crucial to this approach. In particular, it shows that the result found in Theorem 1.1 is superior to the well known results of Kotani theory which prove that $`\gamma (E)>0`$ for almost every $`E`$ in a class of models containing ours. The methods of  also allow to extend Theorem 1.1 to the case where the support of $`\mu `$ is $`\{a,b\}`$, where $`ab`$ are arbitrary real numbers. This requires considerably more technical effort, due to our use of scattering theory at a periodic background, and thereby somewhat obscures the rather simple ideas which we present here in Section 2.

## 2. Positivity of the Lyapunov Exponent

To prove positivity of $`\gamma `$, we need to understand properties of the transfer matrices. The transfer matrix from -1/2 to 1/2 of

$$u^{\prime \prime }+fu=Eu,$$
(9)
is the matrix, $`g(E)`$, for which

$$\left(\begin{array}{c}u(1/2)\\ u^{}(1/2)\end{array}\right)=g(E)\left(\begin{array}{c}u(1/2)\\ u^{}(1/2)\end{array}\right)$$
(10)
for any solution $`u`$ of (9). By $`g_0(E)`$ we denote the corresponding transfer matrix of $`u^{\prime \prime }=Eu`$. Set $`G(E)`$ to be the closed subgroup of $`SL(2,)`$ generated by $`g_0(E)`$ and $`g(E)`$. Let $`P(^2)`$ be the projective space, i.e. the set of the directions in $`^2`$ and $`\overline{v}`$ be the direction of $`v^2\{0\}`$. Note that $`SL(2,)`$ acts on $`P(^2)`$ by $`g\overline{v}=\overline{gv}`$. We say that $`GSL(2,)`$ is strongly irreducible if and only if there is no finite $`G`$-invariant set in $`P(^2)`$.

In order to prove Theorem 1.1 we will consider energies $`E>0`$ and $`E<0`$ separately. Note that $`E=0`$ is contained in $`M`$, so we don’t need to consider it.

The general approach is as follows: We will first prove that $`G(E)`$ is not compact by showing that a sequence of elements has unbounded norm. This argument will be valid for all $`EM^{}`$, where

$`M^{}`$ $`=`$ $`\{(n\pi )^2:n\}\{E=k^2:k>0,b(k)=0\}`$
$`\{\alpha ^2<0:a(i\alpha )a(i\alpha )b(i\alpha )b(i\alpha )=0\}.`$
Once non-compact, then the group $`G`$ is known to be strongly irreducible if and only if for each $`\overline{v}P(^2)`$,

$$\mathrm{\#}\{g\overline{v}:gG\}3,$$
(11)
see . In order to prove that this condition is also satisfied we will in addition have to exclude $`k`$’s which are odd multiples of $`\pi /2`$. We then use Fürstenberg’s Theorem which states, in our context, that if the group $`G(E)`$ is not compact and strongly irreducible then $`\gamma (E)>0`$, see also .

We start with positive energies, i.e. $`E=k^2`$, $`k>0`$, and express the transfer matrices over $`[1/2,1/2]`$ in terms of the Jost solutions (5), at $`\pm k`$. We have

$$g(E)=\left(\begin{array}{cc}u_N(1/2,k)\hfill & u_D(1/2,k)\hfill \\ u_N^{}(1/2,k)\hfill & u_D^{}(1/2,k)\hfill \end{array}\right)$$

where $`u_N`$ and $`u_D`$ are the solutions of (4) satisfying $`u_N(1/2,k)=u_D^{}(1/2,k)`$ $`=1`$, $`u_N^{}(1/2,k)=u_D(1/2,k)=0`$. Writing $`u_N`$ and $`u_D`$ as linear combinations of $`u_+(x,k)`$ and $`u_{}(x,k)=\overline{u_+(x,k)}`$ and setting $`z_\pm (k)=a(k)e^{ik}\pm b(k)e^{ik}`$, we see that

$$g(k^2)=\left(\begin{array}{cc}\text{Re}[z_+(k)]& \frac{1}{k}\text{Im}[z_+(k)]\\ k\text{Im}[z_{}(k)]& \text{Re}[z_{}(k)]\end{array}\right).$$
(12)
Clearly,

$$g_0(k^2)=\left(\begin{array}{cc}\mathrm{cos}k& \frac{1}{k}\mathrm{sin}k\\ k\mathrm{sin}k& \mathrm{cos}k\end{array}\right).$$
(13)
Note for $`k>0`$,

$$\left(\begin{array}{cc}1& 0\\ 0& \frac{1}{k}\end{array}\right)\left(\begin{array}{cc}a& \frac{1}{k}b\\ kc& d\end{array}\right)\left(\begin{array}{cc}1& 0\\ 0& k\end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right).$$
(14)
Thus proving non-compactness for $`G(k^2)`$ is equivalent to proving non-compactness for the group $`\stackrel{~}{G}(k^2)`$ conjugate to $`G(k^2)`$ via (14).

If $`b(k)=0`$ we have that $`|a(k)|=1`$, and hence, we may write $`a(k)=e^{i\varphi }`$. By $`(\text{12})`$, one has then that

$$\stackrel{~}{g}(k^2)=\left(\begin{array}{cc}\mathrm{cos}(k+\varphi )& \mathrm{sin}(k+\varphi )\\ \mathrm{sin}(k+\varphi )& \mathrm{cos}(k+\varphi )\end{array}\right),$$

which is a rotation. The same holds for $`g_0(k^2)`$ by (13). Thus $`\stackrel{~}{G}(k^2)`$ is a group of rotations and thereby compact. This not only excludes an application of Fürstenberg, but shows $`\gamma (k^2)=0`$ as a direct consequence of the definition (3).

Now assume $`b(k)0`$. Set $`a(k)=Ae^{i\alpha }`$, $`b(k)=Be^{i\beta }`$, and see that (6) implies $`A^2B^2=1`$. It follows then that $`A>B>0`$, $`A+B>1`$, and $`AB<1`$. With this notation, one may recalculate

$$\stackrel{~}{g}(k^2)=A\left(\begin{array}{cc}\mathrm{cos}(\phi )& \mathrm{sin}(\phi )\\ \mathrm{sin}(\phi )& \mathrm{cos}(\phi )\end{array}\right)+B\left(\begin{array}{cc}\mathrm{cos}(\psi )& \mathrm{sin}(\psi )\\ \mathrm{sin}(\psi )& \mathrm{cos}(\psi )\end{array}\right),$$

where we have set $`\phi :=k\alpha `$ and $`\psi :=k\beta `$.

We wish to show that a sequence of elements in $`\stackrel{~}{G}(k^2)`$ has unbounded norm. To do so, consider an arbitrary unit vector

$$v(\theta ):=\left(\begin{array}{c}\mathrm{cos}(\theta )\\ \mathrm{sin}(\theta )\end{array}\right).$$

The relation

$$\stackrel{~}{g}(k^2)v(\theta )=Av(\phi +\theta )+Bv(\psi \theta )$$

shows that the image of the unit circle under $`\stackrel{~}{g}`$ is an ellipse, with the choices $`\eta :=\frac{1}{2}(\phi +\psi )`$ and $`\eta ^{}:=\frac{\pi }{2}\frac{1}{2}(\phi +\psi )`$ identifying the semi-major and semi-minor axes as follows:

$$\begin{array}{ccc}\stackrel{~}{g}(k^2)v(\eta )=(A+B)v(\frac{\phi \psi }{2})& \text{and}& \stackrel{~}{g}(k^2)v(\eta ^{})=(AB)v(\frac{\pi }{2}+\frac{\phi \psi }{2})\end{array}.$$

Defining $`R(\theta ):=\stackrel{~}{g}(k^2)v(\theta )`$, a short calculation, using $`A^2B^2=1`$, yields

$$R(\theta )1=2B\left[B+A\mathrm{cos}(2\theta +\phi +\psi )\right].$$

As a consequence, we see that $`R(\theta )1`$ is $`\pi `$-periodic and has exactly two roots in $`[0,\pi )`$. In particular, $`R(\theta )=1`$ if and only if $`\mathrm{cos}(2\theta +\phi +\psi )=\frac{B}{A}`$, which shows that the distance between the zeros of $`R(\theta )1`$ is not equal to $`\frac{\pi }{2}`$: recall $`B0`$ (and $`BA`$). Similarly, $`R(\theta )>1`$ if and only if $`\mathrm{cos}(2\theta +\phi +\psi )>\frac{B}{A}`$ and hence $`|\{\theta [0,\pi ):R(\theta )>1\}|>\frac{\pi }{2}`$. As a result, there exists a compact interval $`K`$ (not necessarily in $`[0,\pi )`$) with $`|K|>\frac{\pi }{2}`$ and $`\stackrel{~}{g}(k^2)v(\theta )>c>1`$ for all $`\theta K`$.

Now, applying $`\stackrel{~}{g}`$ once to $`v(\eta )`$, $`\eta `$ as above, produces a new vector with norm greater than one, but the direction, initially $`\eta `$, is possibly altered. A vector with this new direction may not increase in norm by directly applying $`\stackrel{~}{g}`$ again. However, as long as $`k`$ is not an integer multiple of $`\pi `$, then finitely many applications of $`\stackrel{~}{g}_0(k^2)`$, which is rotation by $`k`$, produces a vector with direction in $`K`$. Once in $`K`$, a direct application of $`\stackrel{~}{g}`$ does increase the norm size uniformly by $`c>1`$, as indicated above. In this manner, we can produce a sequence of elements in $`\stackrel{~}{G}(k^2)`$ with unbounded norm. Thus $`\stackrel{~}{G}(k^2)`$ is non-compact if $`b(k)0`$ and $`k`$ is not an integer multiple of $`\pi `$.

To complete the proof of Theorem 1.1 for positive energies, we note that if in addition $`k`$ is not an integer multiple of $`\pi /2`$, then the free transfer matrix produces three distinct elements in projective space, i.e. (11) is satisfied and Fürstenberg’s Theorem implies $`\gamma (k^2)>0`$.

We finally turn to energies $`E<0`$. In this case $`G(E)`$ is non-compact for every $`E<0`$, since the free transfer matrix

$$g_0(E)=\left(\begin{array}{cc}\mathrm{cosh}(\alpha )& \frac{1}{\alpha }\mathrm{sinh}(\alpha )\\ \alpha \mathrm{sinh}(\alpha )& \mathrm{cosh}(\alpha )\end{array}\right)$$

has unbounded powers, where $`\alpha =\sqrt{|E|}`$. It remains to check (11), which we have to do for $`E(\mathrm{},0)\{\alpha ^2:a(i\alpha )a(i\alpha )b(i\alpha )b(i\alpha )=0\}`$.

Let $`E=\alpha ^2`$ such that $`a(\pm i\alpha )0`$ and $`b(\pm i\alpha )0`$. By (5), this is equivalent to

$$g(E)\overline{v}_\pm \{\overline{v}_+,\overline{v}_{}\},$$
(15)
where $`\overline{v}_\pm `$ are the directions of $`v_\pm :=(1,\pm \alpha )^t`$. If $`\overline{v}\{\overline{v}_+,\overline{v}_{}\}`$, then $`\mathrm{\#}\{g_0(E)^n\overline{v}:n\}=\mathrm{}`$. If, on the other hand, $`\overline{v}\{\overline{v}_+,\overline{v}_{}\}`$, then we use (15) to conclude that an initial application of $`g(E)`$ followed by iteration of $`g_0(E)`$ gives an infinite orbit. This shows (11) and completes the proof of Theorem 1.1.

## 3. Critical energies

In this section we discuss the appearance of critical energies in concrete examples. We thereby illustrate the following points: (i) critical energies with vanishing Lyapunov exponent do indeed exist, (ii) the structure and “size” of the set of energies with $`b(k)=0`$ depends strongly on the concrete example, and (iii) at critical energies at least two different types of non-exponential asymptotics of solutions can be observed.

### 3.1. Zero reflection at degenerate gaps

Our first general observation is that if $`H_{\mathrm{per}}^f`$ is the Schrödinger operator with $`1`$-periodic potential $`V`$ satisfying $`V(x)=f(x)`$, $`1/2x1/2`$ and $`H_{\mathrm{per}}^f`$ has a degenerate gap at energy $`E=k^2`$, i.e. either all solutions of (4) satisfy periodic boundary conditions on $`[1/2,1/2]`$ or all solutions satisfy anti-periodic boundary conditions on $`[1/2,1/2]`$, then $`g(k^2)=I`$ or $`g(k^2)=I`$. By (12) this implies $`(a(k),b(k))=(1,0)`$ or $`(1,0)`$, respectively. Conversely, $`(a(k),b(k))=(\pm 1,0)`$ necessarily requires a degenerate gap. Thus, the possibility of degenerate gaps in periodic potentials is one mechanism which leads to reflectionless energies. Whether all reflectionless energies arise in this way is equivalent to deciding if $`b(k)=0`$ necessarily leads to $`a(k)=\pm 1`$. In Example 2 below we will see that this is not the case.

### 3.2. Example 1

Consider first $`f:=\lambda \chi _{[1/2,1/2]}`$ with $`\lambda `$. For $`E=k^2>\mathrm{max}\{0,\lambda \}`$ we get with $`\alpha :=\sqrt{E\lambda }`$

$`a(k)e^{ik}=`$ $`\mathrm{cos}\alpha +i{\displaystyle \frac{2k^2\lambda }{2k\alpha }}\mathrm{sin}\alpha `$
$`b(k)=`$ $`{\displaystyle \frac{i\lambda }{2k\alpha }}\mathrm{sin}\alpha `$
$`g_0(k^2)=`$ $`\left(\begin{array}{cc}\mathrm{cos}k& \frac{1}{k}\mathrm{sin}k\\ k\mathrm{sin}k& \mathrm{cos}k\end{array}\right)`$
$`g(k^2)=`$ $`\left(\begin{array}{cc}\mathrm{cos}\alpha & \frac{1}{\alpha }\mathrm{sin}\alpha \\ \alpha \mathrm{sin}\alpha & \mathrm{cos}\alpha \end{array}\right)`$
We can now distinguish between the following types of critical energies:

Type 1a: $`k^2`$ such that $`b(k)=0`$: This happens for $`E=k^2=(n\pi )^2+\lambda `$, $`n`$. In this case we have $`g(k^2)=\pm I`$, and it easily follows that the group $`g_0,g`$ is bounded. This means that solutions to (2) are bounded and that $`\gamma (k^2)=0`$.

Type 1b: $`k=n\pi `$, $`n`$: This means $`g_0(k^2)=\pm I`$, and again $`g_0,g`$ bounded, $`\gamma (k^2)=0`$. This type is “dual” to Type 1, since $`k^2=n^2\pi ^2`$ can be viewed as the energies where $`d^2/dx^2+\lambda \lambda \chi _{[1/2,1/2]}`$ is reflectionless with respect to $`d^2/dx^2+\lambda `$.

Type 2: $`E=k^2`$, where $`k=(2n1)\pi /2`$, $`n`$. This leads to a critical energy with a rather different asymptotic behavior of solutions if also $`\alpha =(2m1)\pi /2`$, $`m`$, that is, for specific values of $`\lambda =k^2\alpha ^2=\pi ^2(nm)(n+m1)`$. In this case we have

$$g_0(k^2)=\pm \left(\begin{array}{cc}0& 1/k\\ k& 0\end{array}\right),g(k^2)=\pm \left(\begin{array}{cc}0& 1/\alpha \\ \alpha & 0\end{array}\right).$$

Let $`h_n(\omega ):=g_{2n+1}(\omega )g_{2n}(\omega )`$, where $`g_n(\omega )`$ is the transfer matrix of $`H_\omega `$ from $`n1/2`$ to $`n+1/2`$. Then

$$h_n(\omega )=\{\begin{array}{cc}\pm I\hfill & \text{ with probability }p^2+q^2\text{}\hfill \\ \pm \left(\begin{array}{cc}k/\alpha & 0\\ 0& \alpha /k\end{array}\right)\hfill & \text{ with probability }pq\text{}\hfill \\ \pm \left(\begin{array}{cc}\alpha /k& 0\\ 0& k/\alpha \end{array}\right)\hfill & \text{ with probability }qp\text{}\hfill \end{array}$$

If $`h_{n,1}`$ and $`h_{n,2}`$ are the diagonal entries of $`h_n`$, then the sums of $`\mathrm{log}|h_{n,i}|`$, $`i=1,2`$, give a symmetric random walk (observing that $`(p^2+q^2)\mathrm{log}1+pq\mathrm{log}(k/\alpha )+qp\mathrm{log}(\alpha /k)=0`$). Thus it follows that for every $`\epsilon >0`$ and almost surely

$$\underset{n\mathrm{}}{lim\; sup}n^{1/2\epsilon }\mathrm{log}\underset{k=1}{\overset{n}{}}h_k(\omega )=0$$

and

$$\underset{n\mathrm{}}{lim\; sup}n^{1/2+\epsilon }\mathrm{log}\underset{k=1}{\overset{n}{}}h_k(\omega )=\mathrm{},$$

see e.g. . The interpretation of this is that the hull of the (oscillatory) solutions of (2) asymptotically grows like $`\mathrm{exp}c\sqrt{x}`$ or decays like $`\mathrm{exp}(c\sqrt{x})`$ for some $`c>0`$, where a generalization of the Ruelle-Osceledec Theorem from \[11, Section 8\] is used.

### 3.3. Example 2

Consider $`f:=\lambda (\chi _{[1/2,0]}\chi _{(0,1/2]})`$ with $`\lambda >0`$. For $`E=k^2>\lambda `$, we may again explicitly calculate the transfer matrix

$$g(E)=\left(\begin{array}{cc}g_{11}& g_{12}\\ g_{21}& g_{22}\end{array}\right),$$

and see that,

$`g_{11}=`$ $`\mathrm{cos}(\alpha _{}/2)\mathrm{cos}(\alpha _+/2){\displaystyle \frac{\alpha _{}}{\alpha _+}}\mathrm{sin}(\alpha _{}/2)\mathrm{sin}(\alpha _+/2),`$
$`g_{12}=`$ $`{\displaystyle \frac{1}{\alpha _{}}}\mathrm{sin}(\alpha _{}/2)\mathrm{cos}(\alpha _+/2)+{\displaystyle \frac{1}{\alpha _+}}\mathrm{cos}(\alpha _{}/2)\mathrm{sin}(\alpha _+/2),`$
$`g_{21}=`$ $`\alpha _+\mathrm{cos}(\alpha _{}/2)\mathrm{sin}(\alpha _+/2)\alpha _{}\mathrm{sin}(\alpha _{}/2)\mathrm{cos}(\alpha _+/2),`$
$`g_{22}=`$ $`{\displaystyle \frac{\alpha _+}{\alpha _{}}}\mathrm{sin}(\alpha _{}/2)\mathrm{sin}(\alpha _+/2)+\mathrm{cos}(\alpha _{}/2)\mathrm{cos}(\alpha _+/2),`$
for $`\alpha _\pm :=\sqrt{E\pm \lambda }`$. From this and the boundary conditions of the Jost solution, i.e.

$$g\left(\begin{array}{c}e^{ik/2}\\ ike^{ik/2}\end{array}\right)=\left(\begin{array}{c}a(k)e^{ik/2}+b(k)e^{ik/2}\\ ika(k)e^{ik/2}ikb(k)e^{ik/2}\end{array}\right),$$

one can show that

$$b(k)=0\text{ if and only if }\mathrm{sin}(\alpha _+/2)=0=\mathrm{sin}(\alpha _{}/2).$$

Therefore, if we fix some $`\lambda >0`$, then any $`E=k^2>\lambda `$ for which $`b(k)=0`$ must satisfy

$$\frac{\sqrt{E+\lambda }}{2}=n\pi ,\text{ and }\frac{\sqrt{E\lambda }}{2}=m\pi ,$$

for some $`n,m`$. Thus $`\lambda `$ must have the form

$$\lambda =2\pi ^2(n^2m^2),\text{ with }n,m,$$
(16)
and the corresponding reflectionless $`E`$’s are then

$$E=2\pi ^2(n^2+m^2).$$
(17)
From this we conclude that for most values of $`\lambda `$ (e.g., Lebesgue almost every $`\lambda `$), no reflectionless energy exists above $`\lambda `$. On the other hand, every $`\lambda `$ satisfying (16) has only finitely many reflectionless energies. In fact, if $`N:=\frac{\lambda }{2\pi ^2}`$, then $`N=n^2m^2`$ for only finitely many pairs $`(n,m)^2`$. But this finite number can be arbitrarily large, depending on $`\lambda `$: To find at least $`j`$ exceptional energies, choose

$$\frac{\lambda _j}{2\pi ^2}=N_j=2^{j+1}(2^{j1}+1)$$

and see that for $`0\mathrm{}<j`$ one may choose

$$n_{\mathrm{}}:=\left(2^{\mathrm{}}(2^{j1}+1)+2^{j\mathrm{}1}\right)$$

and

$$m_{\mathrm{}}:=\left(2^{\mathrm{}}(2^{j1}+1)2^{j\mathrm{}1}\right)$$

satisfying $`n_{\mathrm{}}^2m_{\mathrm{}}^2=N_j`$. Proof: Take $`a:=2^{\mathrm{}}(2^{j1}+1)`$ and $`b:=2^{j\mathrm{}1}`$ and see that $`(a+b)^2(ab)^2=4ab`$. Therefore,

$$4ab=2^22^{\mathrm{}}(2^{j1}+1)2^{j\mathrm{}1}=2^{j+1}(2^{j1}+1)=N_j.$$

Recall from our remarks in 3.1 that in some cases, the condition $`b(k)=0`$ implies that $`a(k)=\pm 1`$. To see that this is not always the case, we calculate $`a(k)`$ in this example for such values of $`k`$. For $`\lambda `$ as in (16) and $`E`$ as in (17), one can show that $`a(k)=e^{2\pi i\sqrt{n^2+m^2}}`$ if $`n+m`$ is even and $`a(k)=e^{2\pi i\sqrt{n^2+m^2}}`$ if $`n+m`$ is odd.


format: text

text:
# Magnetically Accreting Isolated Old Neutron Stars

## 1 Introduction

Accretion onto isolated old neutron stars (IONS) through the spherical, Bondi-accretion mode Bondi & Hoyle (1944); Bondi (1952) has been investigated as a means to produce a population of neutron stars detectable in the ROSAT X-ray All Sky Survey (RASS; Treves & Colpi (1991); Blaes & Madau (1993); Madau & Blaes (1994); Zane et al. (1995); Popov et al. (2000)). The estimated number of these sources in the RASS was initially high ($``$$`10^3`$$`10^4`$). However, as observations of higher mean velocities in the radio pulsar population (Lorimer et al., 1997; Hansen & Phinney, 1997; Cordes & Chernoff, 1998, and references therein) were considered, the predicted number of detectable IONS accreting from the ISM decreased dramatically ($``$$`10^2`$$`10^3`$). This is because the Bondi accretion rate is a strong function of the NS velocity through the interstellar medium (ISM)111Here, $`v=\sqrt{V^2+C_s^2}`$, the geometric sum of the spatial velocity of the NS and the sound speed in the accreted medium. The sound speed of the ISM is low ($``$10 km s-1) compared to spatial velocities relevant to the present work ($`v`$$``$100 km s-1), so we assume $`vC_s`$ throughout.

$$\dot{M}_{\mathrm{Bondi}}=4\pi (GM_{\mathrm{NS}})^2\rho v^3$$
(1)
where $`\rho `$ is the local mass density of the ISM. Bondi-accreting IONS are expected to have X-ray luminosities $``$$`10^{31}`$ $`\mathrm{erg}`$$`\mathrm{s}^1`$ and effective temperatures $`kT_{\mathrm{eff}}`$ 200 eV, both dependent upon the accretion rate.

There will be a different accretion mode, however, when the NS magnetic field dominates over the gravitational field in attracting accretion from the ISM. In this mode, the warm ionized medium (WIM) attaches to the NS magnetic field at the magnetosphere, accreting directly onto the NS surface. These magnetically accreting isolated old neutron stars (MAGACs, “magics”) contrast with IONS which accrete through the Bondi-mode, in that the Bondi accretion rate decreases with the NS velocity, while in MAGACs, the accretion rate increases with increasing velocity:

$$\dot{M}_{\mathrm{MAGAC}}=\pi R_\mathrm{M}^2v\rho $$
(2)
where $`R_M=(\mu ^2/(4\pi \rho v^2))^{1/6}`$ is the size-scale for the magnetosphere of a non-rotating NS, where the magnetic energy density from the NS’s magnetic moment ($`\mu `$) is equal to the ram-pressure of the WIM. This gives $`\dot{M}_{\mathrm{MAGAC}}v^{1/3}`$. So, for a given magnetic field strength, Bondi accretion dominates at low velocities, while magnetic accretion dominates at high velocities.

The evolutionary scenario of magnetically accreting NSs has been outlined previously, by Harding & Leventhal (Harding & Leventhal, 1992, HL92 hereafter). HL92 found magnetised NSs first spin down via dipole radiation, followed by a longer spin-down period via the propellor effect. They restricted their attention to sources with surface magnetic fields $``$$`10^{13}`$ G, and found that they are unable to accrete onto the compact object, as the propellor spin-down time is longer than the age of the universe. A similar line of argument is followed in the review by Treves et al. Treves et al. (2000) in the context of Bondi accreting IONS. This limitation does not apply to NSs with magnetic fields $``$$`10^{14}`$ G; a NS with magnetic field $`B=`$$`10^{15}`$ G and velocity $`v`$=300 km s-1, spins down in $``$4$`\times 10^9`$ yr.

On the other hand, this same evolutionary scenario was used to argue for magnetic field decay in the IONS RX J0720.4-3125 Wang (1997); Konenkov & Popov (1997), and for NSs in general Livio et al. (1998); Colpi et al. (1998). We suggest, instead, that high magnetic field IONS have not been discovered due to the fact that strong magnetic field sources accreting from the ISM are spectrally harder (we find below) than Bondi accretors, and have been selected against by observers, or confused with background AGN.

Moreover, since HL92’s work, observations have revealed NSs with considerably stronger magnetic fields than considered by HL92, in the range of $`10^{15}`$ G – the magnetars. Pulse-timing of soft-gamma-ray repeaters Kouveliotou et al. (1998, 1999) confirms their existence as proposed in earlier theoretical work Thompson & Duncan (1993); Duncan & Thompson (1992). Circumstantial evidence suggests that anomalous X-ray pulsars have similarly strong $`B`$-fields Thompson & Duncan (1996); Vasisht & Gotthelf (1997); Heyl & Hernquist (1997); and precision pulsar-timing evidence also implies magnetic fields in excess of $`10^{14}`$ G Kaspi et al. (1999), confirming less-precise multi-epoch timing solutions which nonetheless led to a magnetar interpretation Gotthelf et al. (1999). Static solutions for NSs have been proposed for surface fields up to $``$$`10^{18}`$ G Bocquet et al. (1995); Cardall et al. (2001), although no objects with magnetic fields as high as this have yet been observed. Thus, with theoretical and observational motivation, we reconsider the scenario of HL92, but with particular attention paid to NSs with magnetic fields above $`10^{13}`$ G.

In §2 we describe the model for this population, following the scenario of HL92. In §3, we summarize the properties of MAGACs, discuss observational ramifications and conclude.

## 2 Model

### 2.1 Evolutionary Model

We assume that the neutron star is born with a velocity $`v=v_7`$$`10^7`$ $`\mathrm{cm}\mathrm{s}^1`$ relative to the inter-stellar medium, and a magnetic moment $`\mu =\mu _{33}`$ $`10^{33}`$ G cm3, which corresponds to a dipole field strength at the NS surface (10 km) of $`10^{15}`$ G (appropriate for magnetars; Thompson & Duncan (1993); Duncan & Thompson (1992)). The spin period is initially short (P$`10`$ sec). We assume the magnetic field does not decay. We assume a compact object mass $`M=M_{\mathrm{NS}}1.4`$ $`M_{}`$, radius $`R=R_{\mathrm{NS}}`$$`10^6`$ cm. We assume the NS travels through a spatially uniform proton density $`n_p=0.1n_1`$ cm-3 in the WIM Boulares & Cox (1990); Ferriere (1998a, b). A more detailed examination of the properties of this population, which would include a distribution of NS kick velocities, orbits about a realistic galactic potential, and the galactic gas density and ionization state – such as that performed for Bondi accretors and cooling NSs Popov et al. (2000); Popov (2001) – will be the subject of forthcoming work.

First, HL92 argues that magnetic accretion dominates over Bondi accretion when the magnetospheric radius $`R_\mathrm{M}=4.2\times 10^{12}\mu _{33}^{1/3}v_7^{1/3}n_1^{1/6}\mathrm{cm}`$ is larger than the Bondi accretion radius $`R_{\mathrm{Bondi}}=\frac{2GM_{NS}}{v^2}`$. This inequality results in a velocity limitation for a MAGAC:

$$v_{\mathrm{lim}}>94\mu _{33}^{1/5}n_1^{1/10}M_{\mathrm{NS}}^{3/5}\mathrm{km}\mathrm{s}^1$$
(3)
Electromagnetic dipole radiation pressure inhibits accretion from the ISM until the light cylinder is external to the magnetospheric radius ($`R_{\mathrm{L}.\mathrm{C}.}>R_\mathrm{M}`$), which occurs at a spin period of

$$P_{\mathrm{L}.\mathrm{C}.}=880\mu _{33}^{1/3}v_7^{1/3}n_1^{1/6}\mathrm{s}$$
(4)
The timescale to reach this period through dipole radiation is:

$$\tau _{L.C.}=2.3\times 10^7I_{45}\mu _{33}^{4/3}v_7^{2/3}n_1^{1/3}\mathrm{yr}$$
(5)
where $`I_{45}`$ is the moment of inertia in units of $`10^{45}`$ g cm2. This is a short timescale, compared with the life of the NS. After reaching this period, the light cylinder is at the $`R_\mathrm{M}`$; however, the NS still cannot accrete (or, accretes inefficiently) due to the “propeller effect” Illarionov & Sunyaev (1975) in which matter is spun away by the magnetosphere, which is outside the Keplerian orbit in co-rotation with the neutron star (the co-rotation radius). This matter is ejected from the system, carrying angular momentum with it, spinning down the NS until the co-rotation radius reaches $`R_\mathrm{M}`$, at which point the NS will have a spin period:

$$P_\mathrm{P}=3.9\times 10^6\mu _{33}^{1/2}v_7^{1/2}n_1^{1/4}M_{\mathrm{NS}}^{1/2}\mathrm{s}$$
(6)
This spin period is substantially longer (by three orders of magnitude) than any spin period yet observed from a NS. The propeller spin-down timescale, when the accretion rate is dominated by magnetospheric accretion, is Illarionov & Sunyaev (1975); Davies et al. (1979):

$$\tau _\mathrm{P}=2.3\times 10^{10}I_{45}\mu _{33}^{4/3}v_7^{5/3}n_1^{1/3}\mathrm{yr}$$
(7)
The $`\tau _\mathrm{P}`$ timescale dominates the evolutionary time between NS birth and emergence of the NS as a magnetic accretor. It is comparable to the age of the Galaxy ($``$$`10^{10}`$ yr) velocities already observed among the NS population. but only for high-magnetic field sources ($``$$`10^{15}`$ G). This makes magnetic accretors with the right combination of magnetic field and velocity observable.

Following spin-down via the propeller, matter is accreted onto onto the polar cap of the NS. The maximum accretion rate onto the compact object is (HL92):

$$\dot{M}_{\mathrm{MAGAC}}=\pi R_\mathrm{M}^2v\rho =9.6\times 10^7\mu _{33}^{2/3}v_7^{1/3}n_1^{2/3}\mathrm{g}\mathrm{s}^1$$
(8)
This $`\dot{M}_{\mathrm{MAGAC}}`$ is an upper-limit to the mass accretion rate onto the NS surface, and is potentially much lower (see §2.2). Moreover, the Bondi and MAGAC accretion rates are simply related to the magnetospheric and Bondi radius:

$$\frac{\dot{M}_{\mathrm{MAGAC}}}{\dot{M}_{\mathrm{Bondi}}}=\left(\frac{R_\mathrm{M}}{R_{\mathrm{Bondi}}}\right)^2$$
(9)
If we assume free-fall accretion onto the NS surface, this provides an accretion luminosity of:

$$L_{\mathrm{MAGAC}}=1.8\times 10^{28}\mu _{33}^{2/3}v_7^{1/3}n_1^{2/3}M_{\mathrm{NS}}R_{\mathrm{NS}}^1\mathrm{erg}\mathrm{s}^1.$$
(10)
To estimate the effective temperature, we take the size of the polar cap to be $`R_{\mathrm{cap}}=R_{\mathrm{NS}}(R_{\mathrm{NS}}/R_\mathrm{M})^{1/2}`$. This would seem unusually small – $``$ 50 m for $`B=`$$`10^{15}`$ G. Nonetheless, the effective temperature $`T_{\mathrm{eff}}`$ of this emission is:

$$kT_{\mathrm{eff}}=0.39\mu _{33}^{1/4}n_1^{1/8}M_{\mathrm{NS}}^{1/4}R_{\mathrm{NS}}^1\mathrm{keV}.$$
(11)
Note that $`kT_{\mathrm{eff}}`$ is independent of the NS velocity. The accretion luminosity will remain below the local Eddington luminosity if:

$$\frac{\sigma _{\mathrm{es}}}{m_H}\frac{L_{\mathrm{MAGAC}}}{c\pi R_{\mathrm{cap}}^2}\left(\frac{R_{\mathrm{NS}}}{r}\right)^2<\frac{GM_{\mathrm{NS}}}{r^2}$$
(12)
where $`m_H`$ is the proton mass. The electron scattering cross-section in a strong magnetic field $`\sigma _{\mathrm{es}}`$ is suppressed from the Thomson rate $`\sigma _T`$ by $`\sigma _{\mathrm{es}}/\sigma _T=4\times 10^6(E/1\mathrm{keV})^2(B_{\mathrm{QED}}/B)^2`$, with $`B_{\mathrm{QED}}`$=4.4$`\times 10^{14}`$ G, and photon energy $`E`$ Herold (1979). This result gives the constraint $`\mu _{33}^1n_1^{1/2}<3\times 10^9`$. This will not be violated for NSs with $`B>`$$`10^{13}`$ G unless the ISM is un-physically over-dense ($`n_1`$ $``$$`10^6`$). Thus, the accretion rate will be locally sub-Eddington.

Finally, the ratio of $`\dot{M}_{\mathrm{MAGAC}}`$ to $`\dot{M}_{\mathrm{Bondi}}`$ is:

$$\frac{\dot{M}_{\mathrm{MAGAC}}}{\dot{M}_{\mathrm{Bondi}}}=1.8\mu _{33}^{2/3}v_7^{10/3}n_1^{1/3}M_{\mathrm{NS}}^2$$
(13)
Thus, for the same velocity distribution, MAGACs will be more luminous than Bondi accretors. We discuss the observational ramifications of this in §3.

### 2.2 Comments on the Accretion Luminosity

We have two comments on the accretion luminosity. The first regards the accretion at magnetopause, and the second regards the need for the ISM to be ionized in order to be accreted.

Arons & Lea Arons & Lea (1976) performed detailed analytic description of mass transfer across magnetopause in a $`10^{12}`$ G NS system, fed by a $``$600 km s-1 wind from its companion, but with a substantially higher accretion rate at magnetopause (6$`\times 10^{19}`$ g s-1) than we consider here. Even at this lower magnetic field strength and higher mass accretion rate, they find a cusp forms with a stand-off shock, causing an interchange instability which mediates accretion onto the compact object. While inflow at the cusp increases with increasing density, the interchange instability operates more efficiently and prevents very high densities from arising at magnetopause. Three-dimensional MHD simulations find similar results Toropin et al. (1999).

If the interchange instability of Arons and Lea Arons & Lea (1976) permits accretion across the entire NS surface (rather than just at the magnetic poles), then the luminosity would remain the same, while dramatically decreasing the effective temperature, to:

$$kT_{\mathrm{eff}}=0.0076\mu _{33}^{1/3}v_7^{1/12}n_1^{1/6}M_{\mathrm{NS}}^{1/4}R_{\mathrm{NS}}^{3/4}keV.$$
(14)
Further analysis of accretion onto a magnetic dipole with strength $`10^{13}`$$`10^{18}`$ G is necessary to estimate the magnitude of the effect on the luminosity and spectrum. For now, we adopt $`\dot{M}_{\mathrm{MAGAC}}`$ as an upper-limit.

Secondly, to accrete material via magnetic interaction, the material must be ionized. Throughout, we have therefore assumed that MAGACs accrete from the (ionized) WIM, and not from the neutral ISM. However, one might wonder if the MAGAC’s luminosity is sufficient to ionize the ISM by the time it reaches the magnetosphere (assuming that accretion is already taking place). Ionization of the ISM by an accreting IONS has been examined previously Shvartsman (1971); Blaes et al. (1995). These investigations found that a Bondi-accretion luminosity was capable of ionizing the ISM well beyond the Bondi-accretion radius. However, we have applied a similar approach to MAGACs accreting onto the magnetic polar cap, and find the fraction of a neutral ISM which would be ionized by the MAGAC at $`R_\mathrm{M}`$ to be $`1`$, for a broad range of velocities (0.1$`<v_7<`$10) and densities (0.1$`<n_1<`$100) for magnetic field strengths considered here. Therefore, MAGACs must accrete from the WIM, which has a typical density of $``$0.13 cm-3 for a filling factor of 20% in the plane and a scale height of 1 kpc Ferriere (1998b). This justifies our assumption that MAGACs must accrete from the WIM, and not from the neutral ISM.

### 2.3 Limits on the Observability of MAGACs

Here, we list those effects which will limit the observability of MAGACs, and which determine the luminosity, spectral hardness, velocity, magnetic field, and pulse-period parameter space in which they may be observable. In Fig. 1, we divide the NS $`B`$-field dipole strength, velocity parameter space, on the basis of these limits:

1. $`\dot{M}_{\mathrm{MAGAC}}`$$`>`$$`\dot{M}_{\mathrm{Bondi}}`$. As we find above, this results in a velocity limit (Eq. 3), above which the NS is a MAGAC, and below which the NS is a Bondi accretor (broken line).
2. The propeller phase must be shorter than the age of the galaxy. We use as a limit (Eq. 7) that $`\tau _p<`$1$`\times 10^{10}`$ yr (dotted line), and also include the line for $`\tau _p<`$1$`\times 10^9`$ yr. This age limit also effectively sets a magnetic field strength limit: NSs with $`B<10^{13}`$ G will not spin down in a Hubble time. Thus, if magnetic fields in NSs decay to $`<`$$`10^{13}`$ G in less than a Hubble time, they will not be observed as MAGACs.
3. The pulsar must not escape the Galaxy. We adopt an escape speed limit of 700 km s-1 (cf. Miyamoto & Nagai (1975); Paczynski (1990)), although the exact value depends on where in the galaxy the NS is formed. This limit is indicated by the hashed area at $`v>`$700 km s-1.
4. The NS magnetic field should not exceed the stability limit. We assume a surface dipole field of $`B<`$$`10^{18}`$ G (Bocquet et al. (1995); Cardall et al. (2001) ; dashed-dot line).

We see in Fig. 1 that acceptable magnetic fields above $`10^{14}`$ G, and acceptable velocities range from $``$25 km s-1 to the escape speed from the galaxy. We include lines of constant $`L_{\mathrm{MAGAC}}`$ (for density = 0.1 proton cm-3; Eq. 10); the luminosities of MAGACs are between $`10^{28}`$ and $`10^{31}`$ $`\mathrm{erg}`$$`\mathrm{s}^1`$(again, assuming free-fall accretion from the magnetosphere onto the NS). On the right axis of the figure, we mark $`kT_{\mathrm{eff}}`$, which is independent of velocity, and only weakly dependent upon the density of the ISM ($`n^{1/8}`$). MAGACs have $`kT_{\mathrm{eff}}`$=0.3–2.5 keV. These are harder than observed from any of the identified IONSs to date Treves et al. (2000), perhaps due to an observational bias.

Finally, we note that for values of $`B=`$$`10^{12}`$ G and $`v`$=400 km s-1 which are not unusual in pulsars, such an object will not accrete via the Bondi-mode from the ISM, as the magnetic field is in the propeller range. In addition the source will not spin-down in $`<`$$`10^{10}`$ yr timescale; so it will wander the Galaxy as a dipole-spun down, propelling NS.

## 3 Discussion and Conclusions

Strong magnetic field ($`B>`$$`10^{14}`$ G) NSs go through three stages of evolution: (1) dipole emission; (2) propeller slowdown; and finally, the MAGAC stage (3) magnetic accretion onto the NS surface. MAGACs have the following properties: surface magnetic fields between $`10^{14}`$ G up through the NS stability limit ($``$$`10^{18}`$ G); spin periods between $`10^6`$$`10^8`$ s; X-ray luminosities $`10^{28}`$$`10^{31}`$ $`\mathrm{erg}`$$`\mathrm{s}^1`$, for a WIM density of 0.1 protons cm-3; and effective temperatures $`kT_{\mathrm{eff}}`$=0.3–2.5 keV. The low X-ray luminosities of these sources make it unlikely that the very long spin periods will be detected as pulsations with present X-ray instrumentation, as prohibitively long integrations with frequent monitoring would be required.

At the median velocity for radio pulsars, MAGACs will be $`\times 130`$ more luminous than Bondi accretors (for $`v_7=3`$, a typical median velocity for pulsars; (Lorimer et al., 1997; Hansen & Phinney, 1997; Cordes & Chernoff, 1998, and references therein)). Although the birth-rate of magnetars has been estimated to be $``$10% of isolated pulsars Kouveliotou & et al. (1994); Kouveliotou et al. (1998), they can be observed to $``$10 times the distance and $`130\times `$ the volume of the Galactic disk of Bondi accretors at the radio pulsar median velocity. Magnetic accretors, therefore, may be more prevalent in flux limited X-ray catalogs, such as the ROSAT All-Sky-Survey/Bright Source Catalog (RASS/BSC ; Voges et al. (1999)) than the lower-magnetic field Bondi accretors. Note, however, that Popov et al. Popov et al. (2000) found that the detected RASS/BSC INSs are more naturally interpreted as soft, young, nearby cooling neutron stars, rather than Bondi accretors; in this interpretation, no NSs accreting from the ISM have yet been discovered. The relative number of detectable MAGACs vs. Bondi accretors will ultimately depend on the particulars of the distribution of magnetic fields at birth in NSs, the kick velocities (and their possible dependence on $`B`$), as well as magnetic field decay, which we have neglected herein.

A brief comparison between observations and the predicted properties of this population reveal no strong candidates for MAGACs among known source classes.

Among radio pulsars, there are no known objects with implied $`B`$ and $`v`$ which would place them in the MAGAC $`B,v`$ parameter space. The strongest $`B`$ field radio pulsar is 5.5$`\times 10^{13}`$ G Camilo et al. (2000), is just below the low-$`B`$ limit; its velocity is unmeasured. At a velocity of $``$ 300 km s-1, this NS will spin-down on a timescale of $`\tau _\mathrm{P}`$2$`10^{11}`$ yr, much longer than the age of the Galaxy. Thus, there are no observed radio pulsars in which the magnetic field is clearly strong enough to evolve into a MAGAC. This is possibly due to the relatively short timescale that such objects are dipole emitters ($`\tau _{\mathrm{L}.\mathrm{C}.}`$$``$$`10^6`$ yr).

The observed spin periods of magnetars are $`P`$1-15 sec, which, for the magnetic fields implied ($`10^{15}`$ G), are still within the range where electromagnetic dipole radiation will inhibit accretion. In addition, their X-ray luminosities are a factor $``$1000 larger than those we expect for MAGACs. These objects therefore are not presently MAGACs; however, depending on their spatial velocities, they may evolve into MAGACs.

The (model dependent) age of the isolated NS RX J185635-3754 ($``$1 My; Walter et al. (1996); Walter (2001)) is shorter than $`\tau _P`$, unless the NS has close to $`B`$2$`10^{18}`$ G. With such a high magnetic field, along with the implied velocity ($``$200 km s-1), the source would be within the MAGAC parameter space. Its X-ray luminosity is $``$2$`\times 10^{31}`$ $`\mathrm{erg}`$$`\mathrm{s}^1`$, which is within a factor of a few of the expected luminosity for a MAGAC with this $`B`$ and $`v`$. However, the effective temperature ($`kT_{\mathrm{eff}}=0.057`$ keV is substantially below that expected from a MAGAC with $`B\times 10^{18}G`$ (2.6 keV; Eq. 11). The spectral hardness seems to rule out this object as a MAGAC for the polar-cap accreting scenario. Since all previously observed INSs have effective temperatures $`<`$ 200 eV Treves et al. (2000), it seems similarly unlikely that any known INS is a MAGAC. It is possible, however, that an interchange instability at the magnetosphere permits accretion over a larger fraction of the NS, which results in a softer spectrum, such as observed here (ref. 14).

Discovering these objects in the RASS/BSC requires a similar approach taken in previous work to identify IONS – identifying an X-ray source with no apparent optical counterpart, down to an X-ray-to-optical flux ratio of $`f_x/f_{\mathrm{opt}}`$$``$1000, with confirmation coming from an optical detection at an appropriate value. However, MAGACs can be spectrally harder than their Bondi accreting counterparts ($`kT_{\mathrm{eff}}13`$ keV; cf. Popov (2001)), and searches for these objects must take this into account. Also, given their strong magnetic fields, these objects may have a cyclotron emission line, which may aid identification through X-ray spectroscopy Nelson et al. (1995). Moreover, because we have assumed no magnetic field decay, the absence of MAGACs in the Galaxy, in combination with a population synthesis with high magnetic field SGRs and AXPs, would support field decay models.

We conclude that MAGACs may present a population which are equally prevalent in flux limited X-ray catalogs as Bondi accreting sources. However, no observed examples of this class are yet known.

###### Acknowledgements.

I thank Lars Bildsten, Omer Blaes, David Chernoff and Chris Thompson for useful discussions. I thank Lars Bildsten and Omer Blaes for comments on the draft of this paper, which improved the manuscript. I am grateful to the participants in the ITP Workshop, “Spin and Magnetism in Young Neutron Stars” for stimulating discussions which brought about this work. This research was supported in part by the National Science Foundation under Grant No. PHY99-07949.


text:
# A note on behaviour at an isotropic singularity

## 1 Introduction

Isotropic singularities were introduced into the literature in  as a tool with which to study Penrose’s Weyl Tensor Hypothesis and in an attempt to unify various studies of the initial cosmological singularity. Since then, they have received a lot of attention in the context of this and other hypotheses; in particular, the Vanishing Weyl Tensor Conjecture (that a cosmology admitting an isotropic singularity with initially vanishing Weyl tensor evolves uniquely to a Robertson-Walker universe) has recently been proven for a class of equations of state . The basic idea of an isotropic singularity is that the singular behaviour may be removed by a conformal rescaling, where the conformal factor vanishes at the singularity. The formal definition runs as follows .

A space-time $`(\widehat{},\widehat{g})`$ is said to admit an isotropic singularity if there exists a space-time $`(,g)`$, a smooth cosmic time function $`t`$ defined on $``$ and a function $`\mathrm{\Omega }(t)`$ such that

* $`\widehat{}`$ is the open submanifold $`t>0`$.
* $`\widehat{g}_{ab}=\mathrm{\Omega }^2(t)g_{ab}`$ on $`\widehat{}`$, with $`g_{ab}`$ $`C^3`$ and non-degenerate on an open neighbourhood of $`t=0`$.
* $`\mathrm{\Omega }C[0,b]C^3(0,b]`$ for some $`b>0`$, $`\mathrm{\Omega }(0)=0`$ and $`\mathrm{\Omega }(t)>0`$ for $`t(0,b]`$.
* The limit

$$\lambda =\underset{t0^+}{lim}\frac{\mathrm{\Omega }\mathrm{\Omega }^{\prime \prime }}{(\mathrm{\Omega }^{})^2}1$$

exists, where the prime indicates differentiation with respect to $`t`$.

The low differentiability requirement on $`\mathrm{\Omega }`$ at $`t=0`$, which is identified with the singularity, is required to allow sufficient generality amongst this class of space-times. (If $`\mathrm{\Omega }C^{\mathrm{}}[0,b]`$ and space-time is filled with a perfect fluid with a barotropic equation of state, then the fluid pressure $`p`$ and density $`\rho `$ must obey $`\rho =3p`$ .) The following results may be derived using the conditions imposed upon $`\mathrm{\Omega }`$ and are essential for what follows :

###### Lemma 1

$`\lambda <1`$.

###### Lemma 2

$`lim_{t0^+}(\mathrm{\Omega }^{})^1\mathrm{\Omega }=0.`$

###### Lemma 3

$`(\mathrm{\Omega }^{})^1\mathrm{\Omega }(1\lambda )t`$ as $`t0^+`$.

Such singularities are usually studied in a cosmological context, and so it is generally assumed that there exists a congruence of time-like curves in $`(\widehat{},\widehat{g})`$ representing the aggregate flow of matter in the space-time. We will assume this congruence to be geodesic. So let $`\widehat{v}^a`$ be the unit tangent to a congruence of time-like geodesics of $`(\widehat{},\widehat{g})`$, with typical member $`\widehat{\gamma }`$. We note that as the definition of an isotropic singularity is independent of the matter content of the space-time in question, this congruence need not have the interpretation mentioned above. So we are studying any time-like geodesic (subject to one further restriction given below) which emanates from an isotropic singularity. We can associate with this congruence a congruence of time-like curves in $`(,g)`$ with unit tangent $`v^a=\mathrm{\Omega }\widehat{v}^a`$. The congruence associated with $`\widehat{v}^a`$ is said to be regular if the conformal image $`v^a`$ is $`C^3`$ on an open neighbourhood of $`t=0`$. So let $`(\widehat{},\widehat{g})`$ admit an isotropic singularity, and a regular time-like geodesic congruence with unit tangent $`\widehat{v}^a`$.

Our aim here is to study carefully the idea of the isotropy of the singularity. To this end, we show that every Jacobi field along $`\widehat{\gamma }`$ is crushed to zero length at the same rate as the singularity is approached. Thus physical objects carried along $`\widehat{\gamma }`$ are crushed to zero volume at the same rate in every direction; the crushing effect is isotropic. The rate at which this crushing takes place depends only on the index $`\lambda `$. This is perhaps not surprising; indeed it verifies on a rigorous level part of what one wants to model for such singularities (see Tod’s comments about blobs of fluid in ), but we shall see that the crushing effect arises not just from the fact that $`\mathrm{\Omega }`$ vanishes at the singularity, but that the regularity condition of the existence of the limit $`\lambda `$ and the assumption that $`t`$ is smooth on $`(,g_{ab})`$ are crucial. As a useful by-product, we obtain the general asymptotic behaviour of the conformal factor $`\mathrm{\Omega }`$ as a function of proper time along $`\widehat{\gamma }`$ in the limit as the singularity is approached. We use the metric and curvature conventions of .

## 2 Decomposition and solution of the geodesic deviation equation.

A Jacobi field along $`\widehat{\gamma }`$ is a vector field $`\widehat{\eta }^a`$ orthogonal to $`\widehat{v}^a`$ which satisfies the geodesic deviation equation

$$\widehat{D}^2\widehat{\eta }^a=\widehat{R}_{}^{a}{}_{bcd}{}^{}\widehat{\eta }^c\widehat{v}^b\widehat{v}^d,$$
(1)
where $`\widehat{R}_{}^{a}{}_{bcd}{}^{}`$ is the Riemann tensor of $`(\widehat{},\widehat{g})`$ and $`\widehat{D}=\widehat{v}^a\widehat{}_a`$ is the covariant derivative along $`\widehat{\gamma }`$. The solutions of (1) form a six-dimensional vector space orthogonal to $`\widehat{v}^a`$. To solve (1), we use the assumptions above to identify the singular and non-singular terms of the right hand side. We obtain the asymptotic behaviour of these terms in the limit $`s0^+`$, where $`s`$ is the proper time along an integral curve of $`\widehat{v}^a`$, whose origin is fixed so that the (past) singularity is approached as $`s0^+`$.

By the standard formulae for conformally related metrics ,

$`\widehat{R}_{}^{a}{}_{bcd}{}^{}\widehat{\eta }^c\widehat{v}^b\widehat{v}^d`$ $`=`$ $`\mathrm{\Omega }^2R_{}^{a}{}_{bcd}{}^{}\widehat{\eta }^cv^bv^d`$ (2)
$`+{\displaystyle \frac{1}{4}}(\delta _c^a\mathrm{\Omega }_{bd}g_{bc}\mathrm{\Omega }_d^a\delta _d^a\mathrm{\Omega }_{bc}+g_{bd}\mathrm{\Omega }_c^a)\widehat{\eta }^cv^bv^d,`$
where

$$\mathrm{\Omega }_{ab}=4\mathrm{\Omega }^3_a_b\mathrm{\Omega }+8\mathrm{\Omega }^4_a\mathrm{\Omega }_b\mathrm{\Omega }2\mathrm{\Omega }^4_c\mathrm{\Omega }^c\mathrm{\Omega }g_{ab}.$$

The assumptions made above allow us to simplify the right hand side of (2) considerably. We note first that as $`\widehat{v}^a`$ is assumed to be geodesic, then the relationship between the acceleration vectors of $`\widehat{v}^a`$ and $`v^a`$ which comes from the conformal transformation

$$\widehat{D}\widehat{v}^a=\mathrm{\Omega }^2Dv^a+\mathrm{\Omega }^3D\mathrm{\Omega }v^a+\mathrm{\Omega }^3^a\mathrm{\Omega }$$

yields

$$^a\mathrm{\Omega }=D\mathrm{\Omega }v^a\mathrm{\Omega }a^a,$$
(3)
where $`D=v_a^a`$ and $`a^b=v^a_av^b`$. Then also

$$_a\mathrm{\Omega }^a\mathrm{\Omega }=(D\mathrm{\Omega })^2+\mathrm{\Omega }^2a_ba^b$$

and

$$\widehat{\eta }_a^a\mathrm{\Omega }=\mathrm{\Omega }a_b\widehat{\eta }^b,$$
(4)
and by definitions, $`\widehat{\eta }_a\widehat{v}^a=\widehat{\eta }_av^a=0`$.

The scalar term in the first bracketed term of (2) can be written

$$\mathrm{\Omega }_{bd}v^bv^d=4\mathrm{\Omega }^1\widehat{D}^2\mathrm{\Omega }+2\mathrm{\Omega }^2(\widehat{D}\mathrm{\Omega })^22\mathrm{\Omega }^2a_ba^b,$$
(5)
where we have used $`\widehat{D}=\mathrm{\Omega }^1D`$ when acting on scalars.

Next, a careful calculation relying heavily on (3) and geodesy of $`\widehat{v}^a`$ shows that

$`\mathrm{\Omega }_{bc}\widehat{\eta }^cv^b`$ $`=`$ $`4\mathrm{\Omega }^2Da_b\widehat{\eta }^b.`$ (6)
Finally, the same conditions show that

$`\mathrm{\Omega }_b^a\widehat{\eta }^b`$ $`=`$ $`4\mathrm{\Omega }^2(\widehat{D}\mathrm{\Omega }^av_b+(\widehat{D}\mathrm{\Omega }v^a+a^a)a_b+^aa_b)\widehat{\eta }^b`$ (7)
$`+2\mathrm{\Omega }^2((\widehat{D}\mathrm{\Omega })^2a_ba^b)\widehat{\eta }^a.`$
Thus from equation (2) and equations (4)-(7) the geodesic deviation equation (1) reads

$$\widehat{D}^2\widehat{\eta }^a+\mathrm{\Lambda }_{}^{a}{}_{b}{}^{}\widehat{\eta }^b=0,$$
(8)
where

$$\mathrm{\Lambda }_{}^{a}{}_{b}{}^{}=\mathrm{\Lambda }_{1}^{}{}_{}{}^{a}{}_{b}{}^{}+\mathrm{\Lambda }_{2}^{}{}_{}{}^{a}{}_{b}{}^{}+\mathrm{\Lambda }_{3}^{}{}_{}{}^{a}{}_{b}{}^{},$$

with

$`\mathrm{\Lambda }_{1}^{}{}_{}{}^{a}{}_{b}{}^{}`$ $`=`$ $`\mathrm{\Omega }^2(R_{}^{a}{}_{cbd}{}^{}v^cv^da^aa_bv^aDa_b^aa_b)=:\mathrm{\Omega }^2U_b^a,`$ (9)
$`\mathrm{\Lambda }_{2}^{}{}_{}{}^{a}{}_{b}{}^{}`$ $`=`$ $`\mathrm{\Omega }^2\widehat{D}\mathrm{\Omega }(^av_b+v^aa_b)=:\mathrm{\Omega }^2\widehat{D}\mathrm{\Omega }W_b^a,`$ (10)
$`\mathrm{\Lambda }_{3}^{}{}_{}{}^{a}{}_{b}{}^{}`$ $`=`$ $`\mathrm{\Omega }^1\widehat{D}^2\mathrm{\Omega }\delta _b^a.`$ (11)
We note an essential point for what follows: the definition of an isotropic singularity and the assumption that $`\widehat{v}^a`$ is regular, i.e. that $`v^a`$ is $`C^3`$, implies that the tensors $`U_b^a,W_b^a`$ are finite at the singularity. In the usual way, we can project (8) onto a parallel propagated tetrad, and so treat $`\widehat{\eta }^a`$ as a Euclidean 3-vector and $`\widehat{D}`$ as an ordinary derivative (see section 4.1 of ).

In order to obtain the asymptotic behaviour of solutions of (8) as $`s0^+`$, we must determine the asymptotic behaviour of $`\mathrm{\Omega }`$ and its $`s`$derivatives ($`\widehat{D}\mathrm{\Omega }`$ etc.) in this limit. As a first step towards this, note that

$$\widehat{D}\mathrm{\Omega }=\mathrm{\Omega }^1D\mathrm{\Omega }=\mathrm{\Omega }^1\mathrm{\Omega }^{}Dt.$$

By the assumption that $`t`$ is a smooth cosmic time function throughout a neighbourhood of $`t=0`$ in $`(,g)`$ and the fact that $`v^a`$ is $`C^3`$, we must have that111It is of course possible that $`Dt|_{t=0}=0`$. However the following analysis still holds in this case, as regardless of whether or not $`\dot{t}_0^1\ddot{t}_0`$ converges or diverges, the dominant term in $`L_{}^{a}{}_{b}{}^{}`$ is always $`L_{3}^{}{}_{}{}^{a}{}_{b}{}^{}`$ as can be seen from equations (18-20) in the text.

$$\dot{t}_0:=\underset{s0^+}{lim}Dt>0.$$

(Recall that $`s0^+`$ corresponds to $`t0^+`$.) Thus

$$\widehat{D}\mathrm{\Omega }\dot{t}_0\frac{\mathrm{\Omega }^{}}{\mathrm{\Omega }}.$$
(12)
Similarly, the limit

$$\ddot{t}_0:=\underset{s0^+}{lim}D^2t$$

must exist and be finite, but may be positive, negative or zero. From these observations, we find that

$$\widehat{D}^2\mathrm{\Omega }\mathrm{\Omega }^2\mathrm{\Omega }^{}\ddot{t}_0+(\mathrm{\Omega }^2\mathrm{\Omega }^{\prime \prime }\mathrm{\Omega }^3\mathrm{\Omega }_{}^{}{}_{}{}^{2})\dot{t}_0^2.$$
(13)
Next, consider the term $`\frac{\mathrm{\Omega }^2/\mathrm{\Omega }^{}}{s}`$. Both numerator and denominator vanish in the limit $`s0`$ (see Lemma 2 above), and so we can try to use l’Hopital’s rule to evaluate the limit of the ratio. This gives a finite and hence legitimate result:

$`\underset{s0^+}{lim}{\displaystyle \frac{\mathrm{\Omega }^2/\mathrm{\Omega }^{}}{s}}`$ $`=`$ $`\underset{s0^+}{lim}\left(2\mathrm{\Omega }{\displaystyle \frac{\mathrm{\Omega }^2\mathrm{\Omega }^{\prime \prime }}{(\mathrm{\Omega }^{})^2}}\right)\widehat{D}t`$ (14)
$`=`$ $`\underset{s0^+}{lim}\left(2{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{\prime \prime }}{(\mathrm{\Omega }^{})^2}}\right)Dt`$
$`=`$ $`\dot{t}_0(2\lambda ).`$
The same reasoning may be applied to the quantity $`\mathrm{ln}\mathrm{\Omega }/\mathrm{ln}s`$:

$`\underset{s0^+}{lim}{\displaystyle \frac{\mathrm{ln}\mathrm{\Omega }}{\mathrm{ln}s}}`$ $`=`$ $`\underset{s0^+}{lim}{\displaystyle \frac{\widehat{D}(\mathrm{ln}\mathrm{\Omega })}{\widehat{D}(\mathrm{ln}s)}}`$ (15)
$`=`$ $`\underset{s0^+}{lim}{\displaystyle \frac{s\mathrm{\Omega }^{}Dt}{\mathrm{\Omega }^2}}`$
$`=`$ $`{\displaystyle \frac{1}{2\lambda }},`$
where we have used (14). This yields

$$\mathrm{\Omega }s^{\frac{1}{2\lambda }}.$$
(16)
Recall that $`\lambda <1`$, so that this power is always positive. We note that the relations (14) and (16) hold for all isotropic singularities with a regular congruence of time-like curves: the assumption that $`\widehat{\gamma }`$ is geodesic was not used in the derivation of these equations.

From equations (12)-(16) we can determine the asymptotic behaviour of the singular terms in (8). It is convenient to rewrite (8) in the form

$$s^2\widehat{D}^2\widehat{\eta }^a+L_{}^{a}{}_{b}{}^{}\widehat{\eta }^b=0,$$
(17)
with the obvious defintion of $`L_{}^{a}{}_{b}{}^{}`$ (and its constiuents $`L_{1}^{}{}_{}{}^{a}{}_{b}{}^{}`$ etc. below). The singular coefficient of $`L_{1}^{}{}_{}{}^{a}{}_{b}{}^{}`$ is

$$s^2\mathrm{\Omega }^2s^{\frac{22\lambda }{2\lambda }}.$$
(18)
We note that

$$0<\frac{22\lambda }{2\lambda }<2.$$

The singular coefficient of $`L_{2}^{}{}_{}{}^{a}{}_{b}{}^{}`$ is

$`s^2\mathrm{\Omega }^2\widehat{D}\mathrm{\Omega }`$ $`=`$ $`s^2\mathrm{\Omega }^3\mathrm{\Omega }^{}Dt`$ (19)
$``$ $`{\displaystyle \frac{1}{2\lambda }}s^{\frac{1\lambda }{2\lambda }}.`$
Notice that this term dominates the coefficient of $`L_{1}^{}{}_{}{}^{a}{}_{b}{}^{}`$. The coefficient of $`L_{3}^{}{}_{}{}^{a}{}_{b}{}^{}`$ satisifes

$`s^2\mathrm{\Omega }^1\widehat{D}^2\mathrm{\Omega }`$ $``$ $`s^2\mathrm{\Omega }^3\mathrm{\Omega }^{}\ddot{t}_0s^2(\mathrm{\Omega }^3\mathrm{\Omega }^{\prime \prime }\mathrm{\Omega }^4(\mathrm{\Omega }^{})^2)\dot{t}_0^2`$ (20)
$`=`$ $`s^2\mathrm{\Omega }^3\mathrm{\Omega }^{}\ddot{t}_0{\displaystyle \frac{s^2(\mathrm{\Omega }^{})^2}{\mathrm{\Omega }^4}}\left({\displaystyle \frac{\mathrm{\Omega }\mathrm{\Omega }^{\prime \prime }}{(\mathrm{\Omega }^{})^2}}1\right)\dot{t}_0^2`$
$``$ $`{\displaystyle \frac{1\lambda }{(2\lambda )^2}}{\displaystyle \frac{\ddot{t}_0}{(2\lambda )\dot{t}_0}}s^{\frac{1\lambda }{2\lambda }}.`$
Thus the dominant term in $`L_{}^{a}{}_{b}{}^{}`$ is $`L_{3}^{}{}_{}{}^{a}{}_{b}{}^{}`$, and we can write

$$L_{}^{a}{}_{b}{}^{}=\frac{1\lambda }{(2\lambda )^2}\delta _b^a+s^{\frac{1\lambda }{2\lambda }}M_{}^{a}{}_{b}{}^{},$$

where $`M_{}^{a}{}_{b}{}^{}(s)`$ is continuous at $`s=0^+`$. The $`s`$dependence of the second term here indicates that $`L_{}^{a}{}_{b}{}^{}`$ has a low degree of differentiability, indicating that the solutions of (17) cannot be written as power series in $`s`$. Nevertheless, the leading order behaviour of a fundamental matrix for (17) can be found by looking for a formal solution whose entries are formal logarithmic sums, that is functions of the form

$$f=\underset{j,k=0}{\overset{\mathrm{}}{}}f_{jk}s^{\mu _j}(\mathrm{log}s)^k,$$

where the $`f_{jk}`$ are formal Laurent series in $`s`$ which vanish for sufficiently large $`j+k`$ and the $`\mu _j`$ are complex numbers . It is easily verified that the leading order term in such a series corresponds to a fundamental matrix for the diagonal equation

$$s^2\widehat{D}^2\widehat{\eta }^a+\frac{1\lambda }{(2\lambda )^2}\widehat{\eta }^a=0.$$

A fundamental matrix for this equation is given by

$$y(s)I_3,$$

where $`I_3`$ is the $`3\times 3`$ identity matrix and $`y(s)`$ satisfies the scalar equation

$$s^2y^{\prime \prime }(s)+\frac{1\lambda }{(2\lambda )^2}y(s)=0.$$
(21)
The general solution of this equation is easily found using the method of Frobenius. There are slight variations (the presence of logarithmic terms in the solution) depending on whether or not $`\lambda =2n/(n+1)`$ for some integer $`n`$; however all that concerns us is that the general solution satisfies

$$y(s)cs^\nu ,$$

where $`c`$ is constant and $`\nu `$ (positive) is the smaller of $`1/(2\lambda )`$ and $`(1\lambda )/(2\lambda )`$. Thus the general solution of (21) satisfies $`y0`$ as $`s0^+`$.

Hence there exists a fundamental matrix $`E`$ for (17) which satisfies

$$E(s)y(s)I_3,\text{as }s0^+$$

and $`y0`$ as $`s0^+`$. Therefore every solution $`\widehat{\eta }^a`$ of (17) shrinks to zero as the singularity is approached; the crushing effect of the singularity is the same in every direction. This reinforces the notion of isotropic behaviour at such singularities. We note further that the rate at which the Jacobi fields are crushed depends only on the index $`\lambda `$ and not on any details of the region of space-time from which the singularity is being approached. Thus as one would expect, there is an accompanying degree of homogeneity allied to the isotropic behaviour at the singularity.

## 3 Unexpected behaviour at a singularity which is nearly isotropic

To further elucidate the role played by the defining conditions of an isotropic singularity, we consider an example of a space-time which has an initial singularity with the same broad features as an isotropic singularity, but which does not qualify as such on a technical basis. We show that the approach to this singularity is not accompanied by the isotropic crushing found above for isotropic singularities.

We consider the spherically symmetric conformally flat space-time (so that $`g_{ab}`$ above is the metric tensor $`\eta _{ab}`$ of Minkowski space-time) with conformal factor

$$\mathrm{\Omega }(r,T)=T\mathrm{exp}\left(\frac{r^2}{T^n}\right),$$
(22)
where $`n>0`$ and $`T,r`$ are respectively a global Cartesian time coordinate and the radius function of Minkowski space-time. We take $`T>0`$; then $`T=0`$ is a past singularity of the space-time. Note that since $`\frac{r^2}{T^n}0`$, $`\mathrm{\Omega }0`$ as $`T0^+`$ for all values of $`r`$. For $`T>0`$, $`r=0`$ is the regular centre of the space-time, and is a time-like geodesic $`\widehat{\gamma }`$ with unit tangent

$$\widehat{v}^a=\mathrm{\Omega }^1\delta _T^a.$$

The equations governing the Jacobi fields along such a geodesic are given in . The two tangential Jacobi fields along the radial geodesic $`r=r_0=`$ constant have norms of the form

$$x(s)=c_1R(s)\frac{ds}{R^2(s)}+c_2R(s),$$
(23)
where $`R`$ is the radius function of the space-time, $`s`$ is proper time along the geodesic and $`c_{1,2}`$ are arbitrary constants giving the 2-parameter general solution. Along the geodesic in question, $`R=r_0\mathrm{\Omega }(r_0,T)`$ and $`ds=\mathrm{\Omega }dT`$. Fixing the origin of $`s`$ so that the singularity $`T=0`$ occurs at $`s=0`$, we see that, modulo an irrelevant multiplicative constant,

$$\underset{s0^+}{lim}x(s)=\underset{T0^+}{lim}\mathrm{\Omega }(r_0,T)\frac{dT}{\mathrm{\Omega }(r_0,T)}.$$

Using l’Hopital’s rule, this gives

$$\underset{s0^+}{lim}x(s)=\underset{T0^+}{lim}\frac{T^{n+1}}{T^n+nr_0^2}=0.$$

Notice that this result applies in the case $`r_0=0`$, and so applies to the geodesic $`\widehat{\gamma }`$.

The norm $`y(s)`$ of the radial Jacobi field is governed by the equation

$$\widehat{D}^2y(s)+Fy(s)=0,$$
(24)
where

$$F(r,T)=\mathrm{\Omega }^4(\mathrm{\Omega }_T^2\mathrm{\Omega }_r^2)\mathrm{\Omega }^3(\mathrm{\Omega }_{TT}\mathrm{\Omega }_{rr}),$$

where the subscripts indicate partial derivatives. Using $`ds=\mathrm{\Omega }dT`$, we have

$$\widehat{D}^2y=\mathrm{\Omega }^2\ddot{y}\mathrm{\Omega }^3\mathrm{\Omega }_T\dot{y},$$

where the overdot represents differentiation with respect to $`T`$. Then a straightforward calculation shows that, along $`\widehat{\gamma }`$, equation (24) may be written

$$T^2\ddot{y}T\dot{y}(2T^{2n}1)y=0.$$

If $`n<2`$, the general solution of this equation vanishes at $`T=0`$. We consider the case $`n=2`$. Then the general solution may be written down explicitly;

$$y(t)=c_1T^{1+\sqrt{2}}+c_2T^{1\sqrt{2}}.$$

Thus the general solution, and hence the general radial Jacobi field, blows up at the singularity $`T=0`$. The same holds for $`n>2`$. So the deformation along this geodesic is not isotropic; infinite stretching occurs along the radial direction, and infinite crushing along the tangential directions. Note however that the characteristic volume along $`\widehat{\gamma }`$ that has norm $`V=|x^2y|`$ satisfies

$$VV_0T^{3\sqrt{2}}0.$$

We investigate which of the conditions for an isotropic singularity this example fails to satisfy. We fix $`n=2`$. Then

$`\eta ^{ab}\mathrm{\Omega }_a\mathrm{\Omega }_b=\mathrm{exp}\left({\displaystyle \frac{2r^2}{T^2}}\right)(1+4{\displaystyle \frac{r^4}{T^4}}).`$ (25)
We note therefore that $`^a\mathrm{\Omega }`$ is timelike for all $`(r,T)`$ with $`T0`$, but is degenerate for $`\mathrm{\Omega }=0`$. Thus $`\mathrm{\Omega }=0`$ does not represent a hypersurface within the family of hypersurfaces $`\mathrm{\Omega }=`$constant. However we know that $`\mathrm{\Omega }=0T=0`$, so the singularity may be identified with the space-like hypersurface $`T=0`$. (This is further verified by showing that each point $`(T=0,r=r_r)`$ of the hypersurface is a past endpoint for both an ingoing and an outgoing radial null geodesic of the space-time.)

According to the definition of , the vanishing of $`^a\mathrm{\Omega }`$ on $`\mathrm{\Omega }=0`$ disqualifies this example from membership of the class of isotropic singularities. Similarly, it is not possible to find an appropriate time function $`Z=Z(\mathrm{\Omega })`$ satisfying $`Z(0)=0`$ and $`^aZ0`$. That is, in the language of , there is no conformal gauge choice which allows this example to be considered an isotropic singularity (see also ). This failure to qualify may be clarified by considering how the example fails to qualify on the basis of the definition above, which does not explicitly include the requirement $`^a\mathrm{\Omega }0`$.

First, we note that taking $`t=\mathrm{\Omega }`$, the function $`t`$ is a cosmic time function on $``$, i.e. $`tC^0()`$ and increases along every future-directed causal curve of $``$222This is the definition of . Note the inequivalence to the definition that $`^at`$ is time-like everywhere. The former definition is more general, allowing (i) less smooth functions; (ii) examples such as $`t=T^3`$, where $`T`$ is a global Euclidean time co-ordinate on Minkowski space-time.. This is clear for the open submanifolds $`t0`$, as equation (25) shows. Let $`s^a`$ be tangent to a future-pointing causal curve $`\sigma `$ of $``$. Then $`s^a_at>0`$ for all $`t0`$, but $`s^a_at|_{t=0}=0`$, so certainly $`t`$ is non-decreasing everywhere along $`\sigma `$. But crossing $`t=0`$ is equivalent to crossing $`T=0`$; this increases from negative to positive along $`\sigma `$, which is future directed. Hence $`t`$ also increases from negative to positive along $`\sigma `$. So $`t`$ does indeed increase, in the strict sense, along every future-directed causal curve of $``$; $`t|_{\sigma (\tau _2)}>t|_{\sigma (\tau _1)}`$ for values $`\tau _2>\tau _1`$ of a parameter $`\tau `$ along $`\sigma `$ which increases into the future. In particular, we have $`Dt=1`$ along the geodesic $`r=0`$, and so $`\dot{t}_0=1`$ for this example.

With $`t=\mathrm{\Omega }`$, conditions (i)-(iv) above are trivially satisfied, with $`\lambda =0`$. Also, the geodesic $`\widehat{\gamma }`$ is regular with geodesic conformal image. Thus it is only the failure of $`t:=\mathrm{\Omega }(r,T)`$ to be a smooth cosmic time function on $`(,\eta _{ab})`$ which precludes this example from satisfying all the hypotheses of an isotropic singularity. In fact $`t`$ fails to be smooth ($`C^{\mathrm{}}`$) only at the point $`(r,T)=(0,0)`$. Therefore this example shows the importance of this assumption for isotropic singularities. In terms of the relationship with the definitions of  and , as pointed out by Tod in , the current definition () assumes implicitly that among the allowed cosmic time functions of the definition above, one - $`\overline{t}`$ say - can be chosen so that $`^a\overline{t}0`$ at the singularity. Assuming smoothness for a particular choice $`t`$ would allow one to find the “correct” $`\overline{t}=\overline{t}(t)`$ (this issue may deserve closer attention). The lack of smoothness of our $`t`$ in this example indicates that no such $`\overline{t}`$ exists; the degeneracy of $`\mathrm{\Omega }`$ at the singularity, i.e. $`^a\mathrm{\Omega }=0`$ is gauge-invariant.

## 4 Conclusions

We have studied some technicalities of the definition of isotropic singularities. As one would expect for such a definition, the behaviour at the singularity is indeed isotropic, in the sense that Jacobi fields carried along a geodesic running into such a singularity are crushed uniformly in every direction. This verifies in a rigorous manner what one intends to model with such a definition; a blob of matter would be crushed uniformly to zero volume as the singularity is approached. We have also found homogeneous behaviour at the singularity; the rate at which the Jacobi fields are crushed depends only on the index $`\lambda `$ of the definition above, and not on the region of space-time from which the singularity is approached. Thus although space-times containing isotropic singularities are usually inhomogeneous , the inhomogeneity is washed out in the approach to the singularity.

We have derived this result for time-like geodesics subject to a certain regularity condition. For non-geodesic motions, the deviation equation becomes more complicated (see equation (4.8) of ). However the generalisation of the result to bounded acceleration curves (for example the fluid flow lines in the examples of ) does not appear to be inconceivably difficult: since in this case equation (3) remains true in the limit $`\mathrm{\Omega }0`$, the asymptotic analysis should follow through.

The example of section 3 above shows how these conclusions depend crucially on a technicality of the definition of , namely that the conformal factor is a function of a cosmic time function which is smooth on the background space-time. In our example, the cosmic time function fails to be even $`C^1`$ throughout the hypersurface $`t=0`$. This manifests itself in the fact that regardless of how one tries to choose the conformal factor $`\mathrm{\Omega }`$ and appropriate time function $`Z`$ (see ) for this particular space-time, one finds the degeneracy $`^aZ=0`$ at the singularity.

## Acknowledgement

I am grateful to an anonymous referee for constructive criticism.


text:
# Search for heavy neutrinos mixing with tau neutrinos

## 1 Introduction

In the Standard Model all fundamental fermions have a right-handed component that transforms as an isosinglet under the $`SU(2)_L`$ gauge group except neutrinos, which are observed only in left-handed form. However, heavy neutrinos which are decoupled from $`W`$ and $`Z`$ bosons and hence are mostly isosinglet (sterile) arise in many models that attempt to unify the presently known interactions into a single gauge scheme, such as Grand Unified Theories or Superstrings inspired models . They are also predicted in models trying to solve the problem of baryo- or leptogenesis in the Universe, in many extended electroweak models, such as left-right symmetric and see-saw models . Their masses are predicted to be within the $`\mathrm{GeV}\mathrm{TeV}`$ range. The existence of a light ($`\mathrm{eV}`$ or $`\mathrm{eV})`$ sterile neutrino is expected in schemes that attempt to solve the presently observed indication from atmospheric, solar and LSND experiments that neutrinos are massive, see e.g.  and references therein. More generally one can also look for an isosinglet neutrino with intermediate mass such as in the $`\mathrm{keV}\mathrm{MeV}`$ range. For instance, such neutrinos with masses in the range 1 - 40 $`\mathrm{keV}`$ were recently considered as a candidate for warm dark matter .

If heavy neutrinos exist, many crucial questions arise. For example, for massive neutrinos the flavour eigenstates ($`\nu _e,\nu _\mu ,\nu _\tau ,\mathrm{}`$) need not coincide with the mass eigenstates ($`\nu _1,\nu _2,\nu _3,\nu _4\mathrm{}`$), but would, in general, be related through a unitary transformation. Such a generalised mixing:

$$\nu _l=\underset{i}{}U_{li}\nu _i;l=e,\mu ,\tau ,\mathrm{},i=1,2,3,4,\mathrm{}$$
(1)
could result in neutrino oscillations when the mass differences are small, and in decays of heavy neutrinos when the mass differences are large. The relevant questions are then, do heavy neutrinos mix with the ordinary neutrinos and, if so, what are the elements of the neutrino mixing matrix $`U_{li}`$?

Since there are no firm predictions for $`m_{\nu _4}`$, experimental searches of $`\nu _4`$ have been performed over a wide range of masses. A heavy neutrino in the low mass region (below a few $`\mathrm{GeV}`$) has been searched for mainly in experiments on the leptonic decay of light mesons and in neutrino experiments  resulting in stringent upper limits on $`|U_{li}|^2`$ down to $`10^7`$, thus constraining the mixing amplitudes $`|U_{e4}|^2`$ and $`|U_{\mu 4}|^2`$ for the electron and muon neutrinos, respectively. Concerning the mixing strength $`|U_{\tau 4}|^2`$ for the MeV mass region studied in this experiment a limit can be derived from some earlier papers , see also . Note that stringent limits on $`|U_{\tau 4}|^2`$ in this mass range have been recently obtained from cosmological and astrophysical considerations . For masses in the 3 GeV to $`200\mathrm{GeV}`$ range the LEP experiments have set limits on $`|U_{\tau 4}|^2`$ varying from $`10^5`$ to $`10^1`$ depending on the $`\nu _4`$ mass .

The motivation and purpose of this work is to search for a neutral heavy lepton $`\nu _4`$ which is dominantly associated with the third generation of light neutrinos, $`\nu _\tau `$, via the mixing term $`|U_{\tau 4}|^2`$. If such a particle exists it might be produced in the decay $`D_s\tau \nu _4`$ at the SPS proton target followed by the decay $`\nu _4\nu _\tau e^+e^{}`$ in the NOMAD detector as is illustrated in Figure 1 (see also Section 3). The experimental signature of these events is clean and they can be selected with small background due to the excellent NOMAD capability for precise measurements of the $`e^+e^{}`$ pair direction. An additional motivation for the present study was the time anomaly observed by the KARMEN1 experiment ( see Section 7). Although in this paper we will assume that $`\nu _4`$ has Dirac mass, the application of our results to the Majorana case is straightforward.

## 2 WANF and NOMAD detector

The CERN West Area Neutrino Facility (WANF) beam line  provides an essentially pure $`\nu _\mu `$ beam for neutrino experiments. It consists of a beryllium target irradiated by 450 GeV protons from the CERN SPS. The secondary hadrons of a given sign are focused with two magnetic elements, the horn and the reflector, located in front of a 290 m long evacuated decay tunnel. Protons that have not interacted in the target, secondary hadrons and muons that do not decay are absorbed by a 400 m thick shielding made of iron and earth. The NOMAD detector is located at 835 m from the neutrino target.

The detector is described in Ref. . It consists of a number of sub-detectors most of which are located inside a 0.4 T dipole magnet with a volume of 7.5$`\times `$3.5$`\times `$3.5 m3: an active target of drift chambers (DC) with a mass of 2.7 tons (mainly carbon), an average density of $`0.1\mathrm{g}/\mathrm{cm}^3`$ and a total thickness of about one radiation length ($`1.0X_0`$) followed by a transition radiation detector (TRD) , a preshower detector (PRS), and an electromagnetic calorimeter (ECAL). The PRS is composed of a plane of horizontal and a plane of vertical proportional tubes preceded by a 9 mm thick lead plate ($`1.6X_0`$). The ECAL consists of 875 lead-glass counters, each about 19 $`X_0`$ deep, arranged in a matrix of 35 rows by 25 columns .

A hadron calorimeter (HCAL) and two muon stations are located just after the magnet coils. A plane of scintillation counters, $`V`$, in front of the magnet was used to veto upstream neutrino interactions and muons incident on the detector. Two planes of scintillation counters $`T_1`$ and $`T_2`$ located before and after the TRD were used for triggering purposes. The trigger for neutrino interactions or decays in the DC target was then $`\overline{V}T_1T_2`$.

The electron identification efficiency in NOMAD is provided mainly by the TRD which has an efficiency of more than 90% for isolated electrons of momentum $`150\mathrm{GeV}/\mathrm{c}`$ for a charged pion rejection factor greater than $`10^3`$ .

## 3 Production and decay

As follows from Eq. (1), any source of $`\nu _\tau `$ will produce all kinematically possible massive eigenstates according to the appropriate mixings. In this search the source of $`\nu _\tau (\overline{\nu }_\tau )`$’s that can potentially generate $`\nu _4`$’s is $`D_s(\overline{D}_s)`$ mesons produced in the reaction $`p+BeD_s+X`$ at the proton target and subsequently decaying leptonically: $`D_s\nu _\tau +\tau `$, $`\tau \nu _\tau +X`$ . Up to a $`\nu _4`$ mass of 190 MeV, the mass difference between the $`D_s`$ and the $`\tau `$, $`\nu _4`$’s can originate both from the $`\nu _\tau `$ produced directly in the $`D_s`$ decay and from the $`\nu _\tau `$ produced indirectly in the subsequent $`\tau `$ decay. For a $`\nu _4`$ mass larger than 190 MeV only $`\nu _\tau `$’s produced indirectly in $`\tau `$ decay can contribute. However these indirect $`\nu _\tau `$’s have a lower acceptance at NOMAD and a harder energy spectrum (Figure 3) resulting in a smaller probability to observe the decay $`\nu _4\nu _\tau e^+e^{}`$ in the NOMAD detector. Therefore this search is limited to $`\nu _4`$ masses smaller than 190 MeV. If $`\nu _4`$ is a long-lived particle, the flux of $`\nu _4`$’s would penetrate the downstream shielding without significant attenuation and would be observed in NOMAD through their $`\nu _4\nu _\tau e^+e^{}`$ decays as illustrated in Figure 1.

For neutrino masses below the $`\pi ^0`$-meson mass, $`m_{\nu _4}m_\pi `$, the dominant heavy neutrino visible decay is $`\nu e^+e^{}`$ with a rate which, for small mixing, is given by

$`\tau _4^1\mathrm{\Gamma }(\nu _4\nu _\tau e^+e^{})=K\left[{\displaystyle \frac{(1+g_L^2+g_R^2)G_F^2m_{\nu _4}^5|U_{\tau 4}|^2}{192\pi ^3}}\right]`$ (2)
where $`g_L=1/2+sin^2\theta _W`$, $`g_R=sin^2\theta _W`$, $`K=1(2)`$ for Dirac(Majorana) particles. The corresponding diagram illustrating the dominant contribution from neutral weak currents to this decay mode is shown in Figure 2b. The branching ratio of the visible decay is given by

$$BR(\nu _4\nu _\tau e^+e^{})=\frac{\mathrm{\Gamma }(\nu _4\nu _\tau e^+e^{})}{\mathrm{\Gamma }_{tot}}0.14$$
(3)
where the total rate $`\mathrm{\Gamma }_{tot}`$ is dominated by the $`\nu _43\nu `$ decay channel.

For $`m_{\nu _4}m_\pi `$ the two-body decay channel $`\nu _4\nu _\tau \pi ^0`$ opens up. This mode is phase space favoured and becomes dominant for $`m_{\nu _4}140\mathrm{MeV}`$ .

In this paper we have studied the decay $`\nu _4\nu _\tau e^+e^{}`$ as a possible manifestation of the presence of $`\nu _4`$’s in the neutrino beam. The occurrence of $`\nu _4\nu _\tau e^+e^{}`$ decays would appear as an excess of isolated $`e^+e^{}`$ pairs in NOMAD above those expected from standard neutrino interactions. The decay $`\nu _4\nu _\tau e^+e^{}`$ cannot be distinguished from the anti-neutrino decay $`\overline{\nu }_4\overline{\nu }_\tau e^+e^{}`$ and the result of this search therefore referes to the sum of these two decays.

The spectra of $`\nu _\tau `$’s produced in the Be target by primary protons were calculated using the approach reported in ref.  (see also ). The contribution of protons not interacting in the Be target and interacting in the SPS beam dump at the end of the decay tunnel was also taken into account. The simulated energy spectra of $`\nu _\tau `$’s pointing to the NOMAD fiducial area are shown in Figure 3. The flux of heavy neutrinos, $`\mathrm{\Phi }(\nu _4)`$, can then be expressed as follows:

$`\mathrm{\Phi }(\nu _4){\displaystyle }d\sigma (p+ND_s(\overline{D}_s)+X)/dE_{\nu _\tau }Br(D_s(\overline{D}_s)\tau +\nu _\tau )`$ (4)
$`\stackrel{~}{\lambda }^{1/2}\stackrel{~}{h}^{1/2}|U_{\tau 4}|^2dE_{\nu _\tau }`$
where $`\sigma (p+ND_s(\overline{D}_s)+X)`$ is the $`D_s`$ meson production cross-section , $`Br(D_s\tau +\nu _\tau )`$ is the $`\tau `$-decay mode branching ratio of the $`D_s`$ , and $`\stackrel{~}{\lambda }^{1/2},\stackrel{~}{h}^{1/2}`$ are the decay phase space and helicity factors, respectively .

Once the $`\nu _4`$ flux was known, the next step was to calculate the $`e^+e^{}`$ spectrum based on the differential $`\nu _4\nu _\tau e^+e^{}`$ decay rate including polarization effects (see Figure 4) . The decay electrons and positrons were tracked through the DC target including bremsstrahlung photons, their conversion and multiple scattering in the target. The details of the NOMAD simulation and reconstruction are described elsewhere .

## 4 Data analysis and selection criteria

The search for $`\nu _4\nu _\tau e^+e^{}`$ described in this paper uses the full data sample collected with the $`\overline{V}T_1T_2`$ trigger  during the years 1996-1998. The data correspond to a total number of protons on target (pots) of $`4.1\times 10^{19}`$. The strategy of the analysis was to identify $`\nu _4\nu _\tau e^+e^{}`$ candidates by reconstructing in the DC isolated low invariant mass $`e^+e^{}`$ pairs that are accompanied by no other activity in the detector. The measured rate of $`e^+e^{}`$ pairs was then compared to that expected from known sources.

The following selection criteria were applied:

* two and only two tracks forming a vertex within the DC fiducial volume of $`2.4\times 2.4\times 3.5\mathrm{m}^3`$ equivalent to a mass of 1.97 tons;
* at least one of the two tracks identified as an electron by the TRD (pion contamination probability $`<10^3`$ );
* any additional track or converted photon  in the event were allowed only if their energies were less than 0.4 GeV or if they could be identified as due to bremsstrahlung photons from one of the two electron candidates;
* no $`\gamma `$’s in the ECAL with energy $`E_\gamma >0.4\mathrm{GeV}`$ (or $`E_\gamma >0.3\mathrm{GeV}`$ for $`\gamma `$’s converted in the PRS), which are incompatible with bremsstrahlung photons from the initial tracks;
* total HCAL energy $`<0.4\mathrm{GeV}`$. This cut serves as an HCAL veto and is confirmed by random trigger events and Monte Carlo (MC) studies ;
* the total energy of the pair must be greater than 4 GeV and its invariant mass $`m_{e^+e^{}}<95\mathrm{MeV}`$ to remove background from pairs of particles other than $`e^+e^{}`$.

Only 207 events passed these criteria. At the next step we used a collinearity variable $`𝒞1cos\mathrm{\Theta }_{\nu e^+e^{}}`$, where $`\mathrm{\Theta }_{\nu e^+e^{}}`$ is the angle between the average neutrino beam direction and the total momentum of the reconstructed $`e^+e^{}`$ pair. A cut on this variable allowed a more effective background suppression. A MC simulation of heavy neutrino decays shows (see Figure 5) that the $`\nu _4\nu _\tau e^+e^{}`$ events have $`𝒞<2\times 10^5`$. This was true for a $`\nu _4`$ mass up to $`190\mathrm{MeV}`$.

In order to avoid biases in the determination of selection criteria, a blind analysis was performed. Events in a signal box defined by $`𝒞<2\times 10^5`$ were excluded from the analysis of the data until the validity of the background estimate in this region was established. This was done by verifying that the MC simulation of standard processes reproduced the data outside the box.

The accuracy of the collinearity determination obtained with MC simulations was checked using a $`\nu _\mu CC`$ data sample with an $`e^+e^{}`$ pair from a photon converted in the DC target at a large ( $`100\mathrm{cm}`$) distance from the primary vertex. Figure 6 shows the $`(1cos\mathrm{\Theta }_{e^+e^{}})`$ distribution of such events in the data and simulation, where $`\mathrm{\Theta }_{e^+e^{}}`$ is the angle between the $`e^+e^{}`$ pair momentum and the line joining the primary vertex to the conversion point. The small difference between the MC and data distributions in Figure 6 would result in an overall efficiency correction of less than 6%. However, in order to conservatively account for possible instrumental effects not present in the MC, the MC efficiency was multiplied by the efficiency of reconstructing $`e^+e^{}`$ pairs with a collinearity variable $`𝒞<2\times 10^5`$ in the data sample of Figure 6 ($`75\%`$). Nonetheless, the two distributions are in reasonable agreement at all energies studied. This validates the resolution in the variable $`𝒞`$ (a few mrad in $`\mathrm{\Theta }_{\nu e^+e^{}}`$) predicted by the MC program.

The reconstruction efficiency for the $`\nu _4\nu _\tau e^+e^{}`$ decay in the NOMAD fiducial volume was calculated from the MC simulation as a function of $`e^+e^{}`$ energy in the range 4 GeV to 50 GeV. The MC simulation was used to correct the data for acceptance losses, experimental resolution and reconstruction efficiencies. Two checks using both experimental data and the MC simulation have been performed in order to verify the reliability of the simulation and to estimate the systematic uncertainties in the $`e^+e^{}`$ pair efficiency reconstruction in the energy range predicted by the simulation.

The first method is to select two samples of reconstructed $`\pi ^0`$’s, one in which the two decay photons reach the ECAL, $`N_{2\gamma }^{\pi ^0}`$, and another in which one of the photons converts in the drift chambers, $`N_{\gamma e^+e^{}}^{\pi ^0}`$. For data and MC events the ratio

$$R_{Data,MC}=\left(\frac{N_{\gamma e^+e^{}}^{\pi ^0}}{N_{2\gamma }^{\pi ^0}}\right)_{Data,MC}$$
(5)
was then formed. The value of the double ratio $`RR=R_{Data}/R_{MC}`$ is then a measure of any differences in $`e^+e^{}`$ reconstruction efficiency between the data and MC. The use of $`e^+e^{}`$ pairs from $`\pi ^0`$ decay enhances the purity of the $`e^+e^{}`$ sample.

The method works well mostly for the low energy region, $`E_{e^+e^{}}10\mathrm{GeV}`$, when the $`2\gamma `$ opening angle is relatively large and the distance between the photons in the ECAL is larger than the ECAL cell size. At higher energies the precision of this method is affected by the statistical uncertainties in the number of $`\pi ^0`$’s reconstructed in the $`2\gamma `$ mode, because the resolution on the $`2\gamma `$ opening angle becomes worse and the $`\pi ^0`$ peak is not well identified anymore.

A similar method allowing a more accurate evaluation of the $`e^+e^{}`$ efficiency correction factor at higher energies is based on the inclusive $`e^+e^{}/\gamma `$ double ratio $`RR`$ defined again as $`RR=R_{Data}/R_{MC}`$ with

$$R_{Data,MC}=\left(\frac{N_{e^+e^{}}^\gamma }{N_\gamma ^\gamma }\right)_{Data,MC}$$
(6)
where $`N_\gamma ^\gamma ,N_{e^+e^{}}^\gamma `$ are the numbers of single isolated photons and $`e^+e^{}`$ pairs in the same data sample of $`\nu _\mu CC`$ events used for the collinearity check.

It was found that the two methods agree quite well in the low energy region and yield a correction factor close to 1. However, in the high energy region the $`e^+e^{}`$ efficiency correction factor varied from $`0.7\pm 0.04`$ to $`0.4\pm 0.03`$ depending on the $`e^+e^{}`$ energy.

## 5 Background

The largest contribution to the background is expected from neutrino interactions yielding a single $`\pi ^0`$ with little hadronic activity in the final state. Neutrino interactions occuring in the coil and iron upstream of the DC fiducial volume may yield an isolated $`e^+e^{}`$ pair if a photon from a $`\pi ^0`$ produced in such an interaction converts in the DC and the accompanying particles are not detected.

Because of the large mass of this upstream material the study of this background would require the simulation of a very large number of events resulting in a prohibitively large amount of computer time. Consequently, only about 10% of the required statistics for $`\nu _\mu CC(NC)`$ inelastic reactions were simulated, while other background components, such as $`\nu _eCC`$, coherent $`\pi ^0`$ production, quasi-elastic reactions and $`\nu _\mu e`$ scattering were simulated with a statistics comparable to the number of events expected from these reactions in the data.

The distribution of the variable $`𝒞`$ for the sum of all the MC samples is shown in Figure 5. The plot covers the region $`𝒞<5\times 10^4`$, which is 25 times larger than the size of the signal box. No $`\nu _\mu CC(NC)`$ event is found in this region. The data outside the box, also shown in Figure 5, are consistent with the MC prediction (19 events observed and $`20\pm 4`$ events predicted). The estimate of background from $`\nu _\mu CC(NC)`$ in the signal box is based on the observation that there are no physical reasons for this background to be other than flat in the region $`𝒞<5\times 10^4`$. Two independent methods were used for the background estimation in the signal region.

The first method is based on the MC. The MC background from the fully simulated reactions was found to be $`2.5\pm 0.8`$ events inside the signal region. For the $`\nu _\mu CC`$ and $`\nu _\mu NC`$ events for which only 10 % of the data statistics was generated, no event was found in the full enlarged region $`𝒞<5\times 10^4`$. Assuming that this background is distributed randomly in this interval and taking into account that the simulated sample corresponds to only 10% of data, we estimate the $`\nu _\mu CC`$ and $`\nu _\mu NC`$ background contribution to be $`0_{0.0}^{+0.4}`$ events inside the signal region.

The second method relies on the data themselves. The agreement between the observed and predicted numbers of events outside the signal box allows to conclude that the number of background events from $`\nu _\mu CC`$ and $`\nu _\mu NC`$ processes is negligible. Thus, by extrapolating the 19 observed events to the signal region with the shape of the fully simulated MC events we obtain a second background estimate of $`2.4\pm 0.9`$.

The final background estimates with the two methods are $`N_{bkg}^{\mathrm{MC}}=2.5_{0.8}^{+0.9}(stat)`$ $`\pm 0.6(syst)`$ events from the MC and $`N_{bkg}^{\mathrm{Data}}=2.4_{0.9}^{+0.9}(stat)\pm 0.7(syst)`$ events from the data, thus providing consistent results. The systematic error includes the uncertainties in the number of pot (5%) and in the coherent $`\pi ^0`$ production cross section (25%). In addition, in the second method we also take into account the systematic errors related to the extrapolation procedure. The total systematic uncertainty was calculated by adding all errors in quadrature. In the following we use the background estimate extracted from the data themselves.

## 6 Results and calculation of limits

Upon opening the signal box we have found one event that passes our selection criteria. This is consistent with the expected background and hence no evidence for isosinglet neutrino decays has been found. We can then determine the $`90\%CL`$ upper limit for the corresponding mixing amplitude $`|U_{\tau 4}|^2`$ from the $`90\%CL`$ upper limit for the expected number of signal events, $`N_{\nu _4\nu _\tau e^+e^{}}^{\mathrm{up}}`$. Using the frequentist approach of ref.  and taking into account the uncertainties in the background estimate  we obtain $`N_{\nu _4\nu _\tau e^+e^{}}^{\mathrm{up}}=2.1`$ events. Since the number of observed events is less than the predicted background, this limit is lower than the corresponding sensitivity, $`S_{\nu _4\nu _\tau e^+e^{}}^{\mathrm{up}}=4.0`$. The sensitivity is defined as the average upper limit obtained in the absence of signal events by an ensemble of experiments with the same expected background . The probability to obtain an upper limit of 2.1 or lower is 29%.

For a given flux $`\mathrm{\Phi }(\nu _4)`$, the expected number of $`\nu _4\nu _\tau e^+e^{}`$ decays occuring within the fiducial length $`L`$ of the NOMAD detector located at a distance $`L^{}`$ from the neutrino target is given by

$`N_{\nu _4\nu _\tau e^+e^{}}={\displaystyle \mathrm{\Phi }(\nu _4)exp(L^{}m_{\nu _4}/p_{\nu _4}\tau _4)[1exp(Lm_{\nu _4}/p_{\nu _4}\tau _4)]}`$ (7)
$`(\mathrm{\Gamma }_{e^+e^{}}/\mathrm{\Gamma }_{tot})\epsilon AdE_{\nu _4}|U_{\tau 4}|^4`$
where $`p_{\nu _4}`$ is the $`\nu _4`$ momentum and $`\tau _4`$ is its lifetime at rest, $`\mathrm{\Gamma }_{e^+e^{}},\mathrm{\Gamma }_{tot}`$ are the partial and total mass dependent $`\nu _4`$-decay widths, respectively, and $`\epsilon `$ is the $`e^+e^{}`$ pair reconstruction efficiency. The acceptance $`A`$ of the NOMAD detector was calculated tracing $`\nu _4`$’s produced in the Be-target or beam dump to the detector taking all relevant momentum and angular distributions into account . As an example for a mass $`m_{\nu _4}=33.9\mathrm{MeV}`$, $`A=1.4\%`$ and $`\epsilon =26\%`$. The flux $`\mathrm{\Phi }(\nu _4)`$ is given by Eq. (4).

The final $`90\%CL`$ upper limit curve in the ($`m_{\nu _4};|U_{\tau 4}|^2`$) plane is shown in Figure 7 together with the Big Bang Nucleosynthesis (BBN) and Supernova SN1987a lower limits obtained in ref. . For the mass range $`m_{\nu _4}m_\pi `$ the shape of the NOMAD limit curve is explained by the contribution from the $`\nu _4\nu _\tau \pi ^0`$ decay to the $`\nu _4`$ decay rate: as the $`\nu _4\nu _\tau \pi ^0`$ decay channel opens up, the branching ratio BR($`\nu _4\nu _\tau e^+e^{}`$) drops rapidly as $`m_{\nu _4}`$ increases.

## 7 The KARMEN time anomaly

Our data can also be used to restrict one of the interpretations of the time anomaly observed by the KARMEN experiment . The anomaly is a bump in the time distribution of $`\nu _e`$ and $`\overline{\nu }_\mu `$ induced events which was expected to be well described by the single exponential from muon decays.111This anomaly, seen in the KARMEN1 data, was not confirmed with the KARMEN2 data . A possible explanation of this effect may be found in ref. . The KARMEN Collaboration interpreted the anomaly as being due to an exotic decay of $`\pi ^+`$-mesons into a muon and a 33.9 MeV new fermion $`X`$ and reported a signal curve for pion branching ratio $`BR(\pi ^+\mu ^++X)\times \mathrm{\Gamma }(Xe^+e^{}\nu )`$ as a function of $`X`$-lifetime with a branching ratio $`BR(\pi ^+\mu ^++X)`$ as small as $`10^{16}`$ . This and other hypotheses explaining the anomaly have been recently extensively investigated both theoretically  and experimentally .

Barger et al.  associated the new particle with a 33.9 MeV isosinglet neutrino dominantly mixing with $`\nu _\tau `$. From our limit curve in Figure 7, we have found using Eqs. (2,3) that for $`m_{\nu _4}=33.9\mathrm{MeV}`$ the lifetime of such a neutrino has to be greater than $`10^2`$ sec while the BBN lower limit of ref.  constrain the $`\nu _4`$ lifetime to be less than 0.1 sec.222The SN1987a limit obtained is stronger than the BBN one, however the supernova results are more model dependent. One may imagine a case when the supernova arguments are not applicable, while the BBN ones still apply . Using the KARMEN signal curve  results in a small window around $`BR(\pi ^+\mu ^++X)10^{12}`$ between our result  and BBN one left untested. Note that the recent PSI result on a search for the 33.9 MeV particle in $`\pi ^+\mu ^++X`$ decay corresponds to $`BR(\pi ^+\mu ^++X)<6.010^{10}`$ at 95% $`CL`$ .

## 8 Conclusion

We found no evidence for the existence of a heavy neutrino, $`\nu _4`$, with mass in the range 10 - 190 MeV which mixes dominantly with the $`\nu _\tau `$ and decays into $`\nu _\tau e^+e^{}`$. For the first time an upper limit on the square of the mixing amplitude, $`|U_{\tau 4}|^2`$, is obtained. Windows between our result and the BBN and SN1987a limits are still of interest for further searches.

Acknowledgements

We gratefully acknowledge the CERN SPS accelerator and beam-line staff for the magnificent performance of the neutrino beam. The experiment was supported by the following funding agencies: Australian Research Council (ARC) and Department of Industry, Science, and Resources (DISR), Australia; Institut National de Physique Nucléaire et Physique des Particules (IN2P3), Commissariat à l’Energie Atomique (CEA), France; Bundesministerium für Bildung und Forschung (BMBF, contract 05 6DO52), Germany; Istituto Nazionale di Fisica Nucleare (INFN), Italy; Russian Foundation for Basic Research, Institute for Nuclear Research of the Russian Academy of Sciences, Russia; Fonds National Suisse de la Recherche Scientifique, Switzerland; Department of Energy, National Science Foundation (grant PHY-9526278), the Sloan and the Cottrell Foundations, USA. F.J.P. Soler is supported by a TMR Fellowship from the European Commission. M.C. Gonzalez-Garcia is supported by the Spanish DGICYT under grants PB98-0693 and PB97-1261, by the Generalitat Valenciana under grant GV99-3-1-01 and by the TMR network grant ERBFMRXCT960090 of the European Union.

We thank Prof. V.Matveev for many interesting discussions and comments. We should like to thank Prof. A.D. Dolgov and Dr. D.V. Semikoz for their interest in our results and for many interesting discussions on the astrophysical and cosmological limits.