text:
# Optimal Paths for Spatially Extended Metastable Systems Driven by Noise

## Abstract

The least action principle is exploited as a simulation tool to find the optimal dynamic path for spatially extended systems driven by a small noise. Applications are presented for thermally activated switching of a spatially-extended bistable system as well as the switching dynamics of magnetic thin films. The issue of nucleation versus propagation is discussed and the scaling for the number of nucleation events as a function of the terminal time and other material parameters is computed.

Thermal and other noises play a very important role in the behavior of many systems. No matter how small the noise is, it always has a non-trivial effect on long enough time intervals. Thermal noise has an especially important role for micro- or nano-devices since the energy barrier between metastable states may come close to $`k_BT`$ in such systems and it may in fact limit the smallest size of such devices. One example is found in magnetic recording industry where the superparamagnetic limit, below which the data retention time becomes too short for commercial purposes, is now considered a serious limiting factor for the maximally achievable storage density cow00 . In order to better control such processes, it is necessary to first obtain a detailed quantitative understanding of the effect of the noises in such systems. For bistable systems that are often used in storage and memory applications, issues such as the mean switching time and the dynamics of the switching path have to be addressed. A major challenge is that most systems of interest are spatially extended and switching cannot be described using only a few degree of freedom. For instance, magnetization reversal in micron-sized magnetic films proceeds by nucleation and domain wall motion instead of coherent rotationkogrke00 . Predictions of reversal rates based stochastic coherent rotation models nee55  (stochastic ordinary differential equations) are off by orders of magnitude bra93 . In these cases, one has to deal with the added complexity of stochastic partial differential equations.

Traditionally the methods of choice for a quantitative understanding of the effect of noise have been the Monte Carlo method or direct simulation of the Langevin equation. For small level of noise, however, these methods become prohibitively expensive. Wile the maximal time step for Monte Carlo or Langevin simulations is still restricted by the deterministic part of the dynamics, the effect of the noise is significant only on time scales which are exponentially large in the inverse of noise amplitude. This time scale is the natural one for the noisy dynamics, and it is unaccessible by Monte Carlo or Langevin simulations. For instance, for the magnetic problem (of interest in magnetic recording), the times scale accessible to Monte Carlo or Langevin is on the order of nanoseconds, and the time scale of interest might be years.

In this paper we introduce a new numerical procedure to study spatially extended systems driven by a small noise. We are interested in the effect of rare events over long periods of time. Our method is based on the theory of large deviations which provides a least action principle for the most probable dynamic path frwe98 . Even tough there has already been works using the least action principle for ordinary differential equations (see e.g. ber98 ), extending this in the context of spatially extended systems present non-trivial numerical difficulties. In the examples we consider, namely thermally activated switching in a system with a bistable potential modeled by the Ginzburg-Landau equations in one and two dimensions and sub-micron-sized magnetic thin films, the numerical computations are non-trivial because switching proceeds by nucleation and domain wall motion which are both very localized in space and time, and require using very fine grids. Also, the diffusion terms present in these equations lead to very bad condition numbers, which results in a loss of accuracy as well as an increase of the numbers of iterations necessary for convergence. We introduce an efficient numerical method to overcome these problems which incorporates amongst other techniques the use of suitable preconditioners, efficient minimizations techniques such as the Quasi-Newton method, etc.

Let us begin by considering the simplest continuum model describing the spatio-temporal states of a bistable system, namely the Ginzburg-Landau equation

$$u_t=\delta u_{xx}\delta ^1V^{}(u),$$
(1)
on the interval $`[0,1]`$ with the boundary conditions $`u(0)=u(1)=0`$. We take $`V`$ to be a standard double-well potential, $`V(u)=\frac{1}{4}(1u^2)^2`$. (1) is expressed in appropriate dimensionless variables in which $`\delta `$ is a small parameter which indicates that the reaction term, $`\delta ^1V^{}(u)`$, is fast while the diffusion is slow. (1) can be considered as the gradient flow associated with the energy

$$E[u]=\frac{1}{2}_0^1\left(\delta u_x^2+2\delta ^1V(u)\right)𝑑x.$$
(2)
The dynamics in (1) has two stable equilibrium states, $`u_+`$ and $`u_{}`$, which minimize the energy (2). When $`\delta `$ is small, $`u_\pm (x)=\pm 1`$ except at a thin boundary layer of width $`\delta `$ at $`x=0,1`$ (see Fig. 1).

Now let us add to (2) a small noise term modeling thermal effects:

$$u_t=\delta u_{xx}\delta ^1V^{}(u)+\sqrt{\epsilon }\eta ,$$
(3)
where $`\epsilon `$ is proportional to the temperature of the system and $`\eta `$ is space-time Gaussian white noise with covariance

$$\eta (x,t)\eta (y,s)=\delta (xy)\delta (ts).$$
(4)
The presence of the noise in (3) destroys the long-time stability of the equilibria $`u_\pm `$. For instance, if the initial state is $`u_+`$, there is a finite probability that the system switches to $`u_{}`$ in any time interval $`[0,T]`$. To quantify the probability of such a switching or, in fact, any other event, one can introduce an action functional, $`S_T[u]`$, such that for small $`\epsilon `$, the probability that the solution $`u`$ of (3) be close to a given path $`\phi `$ on $`[0,T]`$ can be estimated by frwe98 ; fajo82

$$\text{Prob}\{u\phi \}\mathrm{exp}(\epsilon ^1S_T[\phi ]).$$
(5)
The action functional corresponding to (3) is given explicitly by

$$S_T[u]=_0^T_0^1\left(u_t\delta u_{xx}+\delta ^1V^{}(u)\right)^2𝑑x𝑑t.$$
(6)
The large deviation principle in (5) allows us to estimate the probability of various events associated with (3) by constrained minimization of the action functional (6). For instance, the probability $`P_T`$ that the system switches from $`u_+`$ to $`u_{}`$ before time $`T`$ satisfies

$$\underset{\epsilon 0}{lim}\epsilon \mathrm{ln}P_T=\underset{u}{\mathrm{min}}\left\{S_{0,T}[u]\right\},$$
(7)
where the minimization in (7) is constrained by the boundary conditions

$$u|_{x=0}=u|_{x=1}=0,u|_{t=0}=u_+,u|_{t=T}=u_{}.$$
(8)
Furthermore assume that the system switches before time $`T`$. Then the minimizer of (7) gives the optimal switching path between $`u_+`$ and $`u_{}`$ during the time interval $`[0,T]`$, in the sense that the probability that the system switches by another path is exponentially smaller in $`\epsilon `$.

To find the optimal path, we numerically minimize the action functional, subject to the boundary condition in (8). We put down a numerical grid on the space-time domain $`(0,1)\times [0,T]`$. The action functional (6) is discretized using finite difference formulas. Numerical optimization is implemented using a Quasi-Newton method, the BFGS method noc99 . In order to speed up convergence, we used the operator obtained from minimizing the linear part of (6), $`_0^T_0^1(u_t\delta u_{xx})^2𝑑x𝑑t`$, as the preconditioner. Similar procedures were used for the other examples below and numerical details are provided in ren .

Figure 1 shows the sequence of profiles of $`u`$ at different times in $`[0,T]`$ for various values of $`T`$ at a fixed $`\delta `$ rem1 . The switching proceeds by nucleation followed by propagation of domain walls. For large $`T`$, the switching proceeds by propagating a domain wall from one boundary to the other. As $`T`$ becomes smaller, however, the number of nucleation events increase. Such switching scenario with multi-nucleations are relevant at finite temperature ren . In Figure 2, where we display the values the action for the various local minimizers, plotted against $`T`$. Figure 3 shows the space-time plot of $`\left(u_t\delta u_{xx}+\delta ^1V^{}(u)\right)^2`$ during a switching event. This corresponds to the minimal (squared) noise necessary to make a switching.

These results can be understood as follows. Consider the critical points of the energy (2), i.e. the solutions of $`0=\delta u_{xx}\delta ^1V^{}(u)`$ with $`u|_{x=0}=u|_{x=1}=0`$. Besides $`u_+`$ and $`u_{}`$, corresponding to the minimizers of the energy, there are also saddle point configurations, denoted by $`u_S`$, with an increasing number of domain walls. For each saddle point, there is a path joining $`u_+`$ to $`u_{}`$ through that saddle point. Our result shows that, for large $`T`$, the switching path crosses the saddle point configuration with minimum energy, i.e the configuration with a single domain wall (see (a) in Figure 1). As $`T`$ is decreased, one might naively think that the switching will occur via the same path followed at a faster rate. This intuition is wrong. For smaller $`T`$, the optimal switching path crosses a saddle point configuration with increasing energy, i.e. more nucleations and therefore more domain walls, giving rise to the cascade of nucleation events seen in Figure 1. The reason is simply that both nucleation and domain wall motion are noise induced. As $`T`$ decreases, at fixed number of nucleations the speed of propagation of the domain wall must increase in order to achieve complete switching during the allowed time $`T`$. This is energy consuming and, at certain critical values of $`T`$ it becomes more favorable to make an additional nucleation.

The same type of cascade in the number of nucleations is also observed if $`T`$ is kept fixed but $`\delta `$ is decreased, because the cost of domain-wall motion increases as $`\delta `$ decreases, whereas the cost of a nucleation stays the same. In fact, a simple scaling argument gives $`n_{}=C(\delta T)^{1/2}`$, for small $`\delta `$ and $`T`$, where $`n_{}`$ is the number of nucleations in the optimal switching path, and $`C`$ is a numerical constant. Furthermore, one has $`lim_{\epsilon 0}\epsilon \mathrm{ln}P_T=C^{}(\delta T)^{1/2}`$, which gives the envelop of the curves in Fig. 2. The detailed analysis of the cascade process of nucleations will be presented elsewhere.

As an example in two dimensions, we minimize

$$S_T[u]=_0^T_\mathrm{\Omega }\left(u_t\delta \mathrm{\Delta }u+\delta ^1V^{}(u)\right)^2d^2x𝑑t,$$
(9)
where $`\mathrm{\Omega }`$ is the unit square, $`\mathrm{\Omega }=[0,1]\times [0,1]`$. We will present results with two different boundary conditions: (i) $`u=1`$ at $`x=0`$ and $`x=1`$, $`u=1`$ at $`y=0`$ and $`y=1`$, and (ii) $`u=0`$ on the edge of the square. In both cases there are two global minimizers of the energy. One minimizer, $`u_+`$, is close to $`u=1`$ except for the boundary layers at $`y=0`$ and $`y=1`$ (case (i)) or at the edge of square (case (ii)). The other minimizer, $`u_{}`$, is close to $`u=1`$ except for the boundary layers at $`x=0`$ and $`x=1`$ (case (i)) or at the edge of square(case (ii)). In the absence of noise, both minimizers are stable equilibrium states.

In Figure 4, we show the time sequences of the switching process for different values of $`\delta `$ at fixed $`T`$. The overall trend is consistent with what was found in the previous example, namely there are more and more nucleation events as $`\delta ,T0`$.

Finally let us consider thermally activated switching of a magnetic thin film, modeled by the Landau-Lifshitz-Gilbert equation, which after suitable nondimensionalization, reads

$$m_t=f(m):=m\times h_{\text{eff}}\alpha m\times (m\times h_{\text{eff}}).$$
(10)
Here $`m`$ is the magnetization distribution (3d vector) which satisfies $`|m|=1`$, $`h_{\text{eff}}`$ is the effective local field, given in terms of the free energy by $`h_{\text{eff}}=\delta E[m]/\delta m,`$ where

$$E[m]=_\mathrm{\Omega }|m|^2d^3x+_\mathrm{\Omega }\varphi (m)d^3x+_{𝐑^3}|u|^2d^3x.$$

Here the three terms represent respectively energies due to exchange, anisotropy, and stray field. $`\mathrm{\Omega }`$ is the region occupied by the magnetic thin film. The potential $`u`$, defined everywhere in space, solves $`\text{div}\left(u+m\right)=0`$, where $`m`$ is extended as $`0`$ outside the sample.

The first term at the right hand side of (10) describes precession of $`m`$ around $`h_{\text{eff}}`$; the second term, which can be written as $`\alpha m\times (m\times h_{\text{eff}})=\alpha (h_{\text{eff}}(mh_{\text{eff}})m)`$, is a damping term whose strength is measured by the parameter $`\alpha `$. We choose physical parameters corresponding to permalloy, and consider a thin square domain of size $`200nm\times 200nm\times 10nm`$.

The dynamics in (10) has equilibrium states corresponding to the minima of the magnetic energy, $`E[m]`$. Thermal noise effects can be represented by modifying $`h_{\text{eff}}`$ as $`h_{\text{eff}}+\sqrt{\epsilon }\eta `$, where $`\eta `$ is a vector white-noise, and will eventually switch the system from one minimum to any other one. In Fig. 5, (m1) and (m4) show two minima where the magnetization points mostly to the right (white color) or to the left (black color), respectively. For small noise, the optimal path for a switching from (m1) to (m4) during the time interval $`[0,T]`$ is obtained by minimizing the action functional corresponding to (10):

$$S_T[m]=_0^T_\mathrm{\Omega }\left|m_tf(m)\right|^2d^3x𝑑t.$$
(11)
The global minimizer of (11) on a long time interval, $`T=450`$, is presented in Fig. 5. Looking at the energy $`E[m]`$ during the switching one sees that two intermediate minima, (m2) and (m3), and three saddle configurations, (s1), (s2), and (s3), are crossed during the switching. The first step of the switching from (m1) to (m2) through (s1) rotates the interior magnetization by $`90^o`$. The second step from (m2) to (m3) through (s2), which is the most expensive step, switches magnetization at the top and bottom end domains. Finally, the third step from (m3) to (m4) through (s3) is similar to the first one.

In conclusion, we have shown that the least action principle can be turned into an efficient numerical procedure for finding the optimal dynamic path in spatially- extended systems driven by small noise, and we have presented applications to bistable systems modeled by the Ginzburg-Landau equation and the magnetic thin films. Our numerical method circumvents the difficulty of having to compute too many time steps in order to observe the relevant events. It is general and can be applied to a variety of problems in physics, chemistry as well as biology.

Weinan E is partially supported by NSF through a Presidential Faculty Fellowship. Eric Vanden-Eijnden is partially supported by NSF Grant DMS-9510356.


format: text

text:
# References

Black holes in the lab?

Steven B. Giddings111email address: giddings@physics.ucsb.edu. On leave from Department of Physics, University of California, Santa Barbara, CA 93106

Department of Physics and SLAC, Stanford University, Stanford, CA 94305/94309

## Abstract

If TeV-scale gravity describes nature, black holes will be produced in particle accelerators, perhaps even with impressive rates at the Large Hadron Collider. Their decays, largely via the Hawking process, will be spectacular. Black holes also would be produced in cosmic ray collisions with our atmosphere, and their showers may be observable. Such a scenario means the end of our quest to understand the world at shorter distances, but may represent the beginning of the exploration of extra dimensions.

Black holes are perhaps the most profoundly mysterious objects in physics. We have long pondered what happens at the core of a black hole. The answer likely involves radically new physics, including breakdown of space and time, and is still beyond the reach of current approaches to quantum gravity such as string theory. Moreover, Hawking’s discovery of black hole radiance and proposal that black holes violate quantum mechanics has led us to the sharp paradox of black hole information, which drives at the very heart of the problem of reconciling quantum mechanics and gravity. There is no clear way out:222For reviews, see . information loss associated with breakdown of quantum mechanics apparently leads to disastrous violations of energy conservation; information cannot escape a black hole without violating locality; and the third alternative, black hole remnants, lead to catastrophic instabilities. String theorists have recently investigated the second alternative, via holography, but the jury is still out as no one has managed to understand how holographic theories can reproduce the approximately local physics that we see in our everyday world.

Experimental clues to the physics of black hole decay would be welcome. Unfortunately, manufacture of microscopic black holes apparently requires scattering energies above the four-dimensional Planck mass, $`M_410^{19}`$ GeV, placing this possibility far in our future.

However, recently there has been a revolution in thinking about the relationship between the Planck scale and the weak scale, $`M_W1`$TeV. The longstanding “hierarchy problem” is to explain the large ratio of these; one would naturally expect $`M_WM_4`$. The new idea is that the weak scale and the fundamental Planck scale, $`M_P`$, are indeed the same size, but four-dimensional gravity is weak (hence $`M_4`$ is large) due to dilution of gravity in large or warped extra dimensions.

Specifically, for a general Poincaré-invariant metric

$$ds^2=e^{2A(y)}dx^{\mu 2}+g_{mn}(y)dy^mdy^n,$$
(1)
where $`x^\mu `$ are the four dimensions we see, $`y^m`$ are the extra $`n`$ dimensions, and the function $`A(y)`$ is the warp factor, $`M_4`$ is given by

$$\frac{M_4^2}{M_P^2}=M_P^nd^ny\sqrt{g}e^{2A}.$$
(2)
Either large volume or large $`A`$ produce a big ratio between $`M_P`$ and the observed Planck scale. The hierarchy problem morphs into that of explaining why the extra dimensions are large, or highly warped. For $`n=2`$$`6`$ their size ranges from $`mmfm`$. This presents a conflict with precision measurements of the gauge forces which have now reached the scale of $`10^3fm`$, but we are saved by the brane world idea, which follows naturally from string theory: gauge forces, corresponding to open strings, propagate on a brane within the extra dimensions, whereas gravity, which is always transmitted by closed strings, propagates in all of the dimensions. Such scenarios go by the name “TeV-scale gravity.”

If TeV-scale gravity describes nature, the consequences are astounding. We will begin to explore quantum gravity, and possibly string theory, at accelerators in the relatively near future. Indeed, model independent bounds based on present experiments merely indicate $`M_P>800`$GeV; TeV-scale gravity could be the physics of the Large Hadron Collider (LHC).

The most generic and spectacular result of such a scenario would be the production of black holes in particle accelerators. LHC will collide protons, which are aggregates of partons (quarks and gluons). If $`M_P𝒪(TeV)`$, then parton collisions with significantly higher center-of-mass energy $`E`$ should produce black holes; we consider $`E>5M_P`$, such that the Bekenstein-Hawking entropy $`S_{BH}>25`$ (for $`n=6`$) in order to ensure that these are close to having a classical description. In TeV-gravity/brane world scenarios, there are two other important approximations. The first is that the gravitational field of the brane can be neglected, and is valid for black holes heavy compared to $`M_P`$, and the second is to consider black holes small as compared to the scales of the extra dimensions. We may then effectively discuss black holes in $`4+n`$ dimensional flat space, as studied in .

The first question in these scenarios concerns the production rate for black holes. Two ingredients are needed: the parton density in a proton, which is approximately known, and the cross-section for two partons to form a black hole, which is not. This cross-section may, however, be estimated. Arguments along the lines of Thorne’s hoop conjecture indicate that a black hole forms when partons collide at impact parameter $`b`$ that is less than the Schwarzschild radius $`r_h`$ corresponding to $`E`$. This would suggest a parton-parton cross-section of the form

$$\sigma r_h^2(E)E^{\frac{1}{n+1}}.$$
(3)
However, until now the high-energy gravitational collision problem has been little studied. In 1974, Penrose argued that black holes form in zero impact parameter collisions, and this work was extended by D’Eath and Payne, but the problem at non-zero $`b`$ had not been systematically treated. In , Doug Eardley and I recently revisited this problem, and in particular showed that in four dimensions, a trapped surface forms in collisions with impact parameter $`b<1.6E`$, very close to the naïve expectation of $`b<2E`$. Furthermore, in higher dimensions we reduced the problem of finding a trapped surface to a higher-dimensional analog of the Plateau problem, which we expect to have a solution – work on this continues.

Using the estimate (3), one readily finds an impressive result: for $`M_P1`$ TeV, the LHC will produce black holes with masses larger than $`5M_P`$ at the rate of about one per second. This would qualify LHC to be called a black hole factory.

Black holes will then decay leaving spectacular signatures. The first stage of their decay is purely classical, and involves the rather asymmetric initial black hole settling down to a hairless spinning black hole, radiating its multipole moments. We call this stage “balding.” An important open problem is to determine how much energy is left in the black hole at the end of this stage; rough estimates based on the size of the initial trapped surface – which can only grow – and extrapolation of  suggest that this energy is around $`1540\%`$ of the initial energy $`E`$. We hope that improvement of numerical relativity or perturbation methods eventually give us a better characterization of this stage.

Hawking’s calculation then becomes relevant. As in the decay of four-dimensional black holes, one expects the black hole to first shed its spin, radiating particles preferentially in the equatorial plane, in a “spin-down” phase. Extrapolation of Page’s four-dimensional results suggest a mass loss of perhaps $`25\%`$ in spin-down. An important problem is to redo Page’s analysis in the higher-dimensional setting.

Spin-down leaves a Schwarzschild black hole which continues to evaporate through the “Schwarzschild phase.” The instantaneous energy distribution is thermal, and may be integrated to find an overall spectrum. This phase should represent perhaps $`75\%`$ of the black hole decay energy.

When the black hole reaches the Planck size, we confront the profound mystery we began with: what effects govern the final decay, what do they tell us about the nature of quantum gravity, and what happens to information? Exploration of this “Planck phase” is beyond present theoretical technology, which makes the prospect of experimental results all the more tantalizing.

Products of these stages should stand out in accelerators. In particular, a black hole should produce of order $`S_{BH}`$ energetic primary particles – leptons, quarks, gluons, etc. – in its Schwarzschild phase. These will be radiated roughly isotropically, with characteristic spectra and species ratios, as predicted by Hawking’s calculations and numerical extensions. These events should not be masked by backgrounds from any known extrapolation of the standard model – they are very unique.

This scenario has other interesting consequences. First, we know that cosmic rays hit our atmosphere with center-of-mass energies exceeding energies accessible at LHC: if TeV-scale gravity is correct, black hole events have peppered our upper atmosphere throughout earth’s history. If ultrahigh-energy neutrino cosmic ray fluxes are sufficiently strong, these events may even be observable in the next round of cosmic ray observatories.

It also appears likely that, through the AdS/CFT correspondence, relativists might tell particle theorists something about QCD. I have recently argued that black hole formation or other strong gravity effects in anti-de Sitter space are dual to the physics that saturates the Froissart bound for hadron cross-sections

$$\sigma \mathrm{ln}^2E.$$
(4)
Important questions remain regarding the structure and stability of such black holes.

Finally, consequences for the future of high energy physics are profound. Humanity has pursued a longstanding quest to understand physics at shorter and shorter distances. In quantum gravity we don’t know that distances shorter than the Planck scale exist, but this is a question that should be addressed experimentally. However, once we start making black holes this appears impossible. Any attempt at collisions that can probe shorter distances will be cloaked inside an event horizon, and all that will be seen in our detectors is the products of the black hole decay – there is apparently no way to directly observe the short-distance physics taking place inside the black hole. Black hole production represents the end of short-distance physics. However, it doesn’t necessarily spell a dismal future for high-energy experiments. As we produce bigger black holes, they reach off the brane that is our world and offer us a way to explore the geometry and other features of the extra dimensions. High energy physics can become the study of the geography of extra dimensions.

Acknowledgments

I’d like to thank my collaborators D. Eardley, E. Katz, and S. Thomas for the opportunity to explore these fascinating ideas together. This work was supported in part by the Department of Energy under Contract DE-FG-03-91ER40618, and by David and Lucile Packard Foundation Fellowship 2000-13856.


text:
# Microelectromagnets for Trapping and Manipulating Ultracold Atomic Quantum Gases

## Abstract

We describe the production and characterization of microelectromagnets made for trapping and manipulating atomic ensembles. The devices consist of 7 fabricated parallel copper conductors 3 $`\mu `$m thick, 25 mm long, with widths ranging from 3 $`\mu `$m to 30 $`\mu `$m, and are produced by electroplating a sapphire substrate. Maximum current densities in the wires up to $`6.5\times 10^6`$ A cm-2 are achieved in continuous mode operation. The device operates successfully at a base pressure of 10-11 mbar. The microstructures permit the realization of a variety of magnetic field configurations, and hence provide enormous flexibility for controlling the motion and the shape of Bose-Einstein condensates.

The combination of ultracold atoms and miniaturized magnetic traps  suggests fascinating scenarios for integrated atom optic devices including waveguides, interferometers, resonators, and quantum gates . The current progress in reducing the production time of Bose-Einstein condensates below 1 s  and the recently achieved significant simplification of the magnetic apparatus  adds to the vision of a future quantum technology with coherent matter waves analogous to today’s highly developed laser technology. One of the crucial prerequisites is the preparation of atomic clouds in steep potentials. As a particularly attractive approach, magnetic fields generated by miniaturized current conductors can be used . In addition, the trapping potentials can be spatially structured on a micron length scale and temporal modulations are possible on the scale of microseconds. By using miniaturized conductors, the trapping potential is located close to the surface of the substrate. At a distance of about 100 $`\mu `$m from the surface loss mechanisms, heating and decoherence effects within the cloud are predicted to become significant . Earlier experiments succeeded in loading atoms into microtrap potentials that were generated by current conductors on the scale of several 10 $`\mu `$m . Microtraps based on smaller conductors have been demonstrated electrically , and single atoms have been successfully guided in a micro magnetic channel . In this Letter we describe the realization of magnetic microtraps with ultra thin conductors that reach the limit of surface-atom-interaction. Details of the fabrication procedure, and the chemical and electrical properties are given. The particular microstructure described has been used for generating highly anisotropic Bose-Einstein condensates . We give a brief summary of the experiments and conclude with an estimation of the field configurations which can be reached with this device.

The microelectromagnet consists of a parallel configuration of copper wires with widths ranging from 3 $`\mu `$m to 30 $`\mu `$m. The thickness of the galvanically deposited copper was of the order of a few microns, while the length of the total structure was approximately 25 mm.

The fabrication of thin metallic wires with such high aspect ratios and excellent conducting properties requires the use of sophisticated microfabrication techniques. Initially, a high-quality optical mask was prepared using electron-beam lithography, which was then employed to transfer the wire pattern onto the smooth (residual roughness $`1\mu `$m) sapphire substrate. In order to guarantee good adhesion and to provide a contact layer for the subsequent galvanic processing the substrate was also coated with a thin (7 nm Cr followed by 120 nm Cu) metallic layer by thermal evaporation. The wire fabrication itself was performed in a conventional acidic CuSO4 galvanic bath with an additional organic brightener which produces copper layers of high-purity ($`99.99\%`$) and a density of 8.9 g cm-3. Best results were obtained using the technique of pulsed plating where the current was modulated with a square wave of frequency 1 kHz and current values of 360 mA and -40 mA for the positive and negavtive excursions respectively. This technique is well suited for the fabrication of smooth metallic layers as evidenced by the finished structure shown in Fig. 1.

For the chosen cathode area (ca. 44 cm2) the above conditions correspond to a current density of approximately 0.8 A dm-2 (and -0.09 A dm-2 for the negative half-cycle), and to an electroplating rate of 1 nm s-1. After removal of the photoresist the metallic contact layer was chemically etched with nitric acid for copper removal followed by potassium ferricyanide for the chrome removal. The quality of the electroplated layer and substrate was confirmed via EDX measurements which show neither the presence of magnetic impurities nor significant contamination.

A test structure with a 1.6 $`\mu `$m copper layer was used to determine the critical current at which the conductors are thermally damaged. The measurements were carried out at room temperature in air with the microstructure mounted on a copper heat sink. The electrical contact to the microstructure is accomplished by bonding 3 aluminium wires with 20 $`\mu `$m diameter and a length of about 700 $`\mu `$m onto each contact. Roughly 10 s after changing the current in the conductors a steady state temperature is reached, and the voltage drop along the wires can be recorded. Thus the current was increased every 10 s in steps of 10 mA for the 3 $`\mu `$m-conductor and in steps of 50 mA for 11 $`\mu `$m- and 30 $`\mu `$m-conductors. The maximum stable current is marked by the last data point in Fig. 2, which shows the device resistance as a function of the current.

The corresponding maximum current densities are 3.1$`\times 10^6`$ A cm-2, 4.8$`\times 10^6`$ A cm-2, and 6.5$`\times 10^6`$ A cm-2 for the 30 $`\mu `$m, 11 $`\mu `$m, and 3 $`\mu `$m conductors respectively. The maximum current density increases with decreasing width, presumably because the geometry of thinner conductors favours a more efficient heat contact with the substrate. The resistance increases by a factor of 2.5 and of 2.3 respectively which corresponds to a temperature increase of about 350 K, determined using the thermal coefficient of resistance for bulk copper ($`\alpha =0.004`$K-1). We therefore estimate that the average temperature of the wires at breakdown is roughly 650 K. Near imperfections, which may be caused by unintentional variations of the conductor cross section or by reduced adhesion to the substrate, we expect the temperature to exceed the average value and to lead to breakdown when the local temperature reaches the melting point of bulk copper (1350 K). The temperature behaviour near the maximum current is similar for the 11 $`\mu `$m and 30 $`\mu `$m conductors. The resistance of the 3 $`\mu `$m conductor increases only by a factor of 1.5 at the maximum current and the average temperature is increased to only 450 K. This can be explained by a stronger influence of the imperfections. The maximum current densities measured in our test experiments at room temperature are similar to the critical current densities in high-Tc superconductors at liquid nitrogen temperature  and fall short only by a factor 15 of the values obtained by Drndic et al.  with gold conductors on sapphire substrates at liquid helium temperature. The microelectromagnets are used under ultra-high vacuum conditions and can be cooled with liquid nitrogen. A significant resistance reduction and consequent increase in critical current density is observed. In addition, the cooling reduces outgassing of the surface and results in longer lifetimes of the trapped atoms. Nevertheless the longest lifetimes observed in recent experiments  depend critically upon the separation between substrate and condensate; for separations of 200 $`\mu `$m and 20 $`\mu `$m lifetimes of 15 s and 700 ms respectively were found.

The configuration of seven parallel conductors with different widths combined with the extreme length of the microstucture of 25 mm permits the definition of a variety of magnetic trapping potentials . Besides the usual magnetic traps used for Bose-Einstein condensation, very elongated waveguides can also be obtained by changing the current in the microstructure. Furthermore, a parallel set of waveguides can be formed and may be merged and separated by varying the strength of the offset field perpendicular to the waveguides . Such configurations permit the realization of on-chip interferometers, and are currently the subject of investigation.

In order to achieve the Bose-Einstein condensation  one conductor is driven for example with an current source of 0.2 A, and an external bias field of 4 G is applied. The resulting gradient in this instance amounts to 400 G cm-1 and the center of the trap is located 100 $`\mu `$m above the surface of the microstructure. With a superimposed offset field of 1 G aligned parallel to the conductor the curvature of the magnetic field magnitude in the transverse direction amounts to 16 kG cm-2 corresponding to a radial (transverse) oscillation frequency of $`\omega _r=2\pi \times 510`$ s-1 for 87Rb atoms in the $`|`$F=2,mF=2$`>`$ ground state. The axial (longitudinal) confinement of the atoms is separately controlled by an additional Ioffe-type electromagnet , which is adjusted such that the axial oscillation frequency amounts to $`\omega _a=2\pi \times 14`$ s-1. The atoms are loaded into the microtrap by adiabatic compression  and subsequently cooled below the critical temperature by forced rf-evaporation. Condensation is reached with 106 atoms at a temperature of 1 $`\mu `$K. This standard trap configuration can be continuously changed into a much steeper potential with an elongated waveguide of large aspect ratio. In the particular experiment shown in Fig. 3 the axial oscillation frequency of the trap was reduced by a factor of 10 on a timescale of 400 ms. During the last 100 ms the trap center was shifted by 1 mm such that the condensate is free to propagate along the waveguide. This motion is clearly visible, and observations over longer time scales show oscillations of the center of mass of the condensate at the trap frequency of $`2\pi \times 1.4`$ s-1. The distance to the surface and the radial oscillation frequency were tuned to 270 $`\mu `$m and $`2\pi \times 100`$ s-1 accordingly.

In the most extreme potential configuration obtainable with our microstructure, defined by the maximum current density and the size of the thinnest conductor, the magnetic field gradient is expected to reach 500 kG cm-1. With a current of 200 mA in the center conductor and -325 mA in both of the two neighboring conductors the trap is located 3 $`\mu `$m above the surface and has a depth of 133 G corresponding to a temperature of 8 mK for 87Rb atoms. With an axial offset field of 1 G, the radial curvature amounts $`2.5\times 10^{11}`$ G cm-2 corresponding to an oscillation frequency of $`\omega _r=2\pi \times 600,000`$ s-1. By choosing a small axial oscillation frequency (for example $`\omega _a<2\pi \times 5`$ s-1), an aspect ratio of $`\omega _r/\omega _a>`$ 100,000 can be achieved. To date, the properties of cold atomic ensembles have not been investigated under such extreme conditions due to the lack of an appropriate device. In our most recent experiments we observe periodic substructures in the spatial distribution of an elongated condensate at a surface separation below 150 $`\mu `$m, which is reported in detail elsewhere . The obvious causes for such modulated magnetic or electrostatic potentials have been excluded, and we speculate that their origin lies in the magnetic surface properties of the conductors.

In conclusion, we have demonstrated a device for the manipulation of Bose-Einstein condensates based on microfabricated electromagnets. The microfabrication technique employed allows the construction of even more sophisticated micromagnets, and opens new possibilities for integrated atom-optical devices with ultracold atoms. Novel sensors for rotation and for precision measurement of the gavitational field may become feasible.

###### Acknowledgements.

This work is supported in part by the Deutsche Forschungsgemeinschaft under grant no. Zi 419/3-1. We thank S. Raible and R. Frank for technical assistance, and HighFinesse for providing ultra stable current sources for driving the microelectromagnet.


text:
# MIKA: a multigrid-based program package for electronic structure calculations

## I Introduction

The goal of computational materials science and also that of modeling of nanoscale man-made structures is to calculate from first principles the various chemical and/or physical properties. This requires the solution of the electronic (and ionic) structures of the system in question. The density-functional theory (DFT)  makes a huge step towards this goal by casting the untractable problem of many interacting electrons to that of noninteracting particles under the influence of an effective potential. However, in order to apply DFT in practice one has to resort to approximations for electron exchange and correlation such as the local-density approximation (LDA) or the generalized-gradient approximation (GGA). Moreover, in the case of systems consisting of hundreds or more atoms it is still a challenge to solve numerically efficiently for the ensuing Kohn-Sham equations.

We have developed a real-space multigrid method called MIKA (Multigrid Instead of the K-spAce) for the numerical solution of the Kohn-Sham equations . In real-space methods, the values of the wave-functions and potentials are presented using three-dimensional point grids, and the partial differential equations are discretized using finite differences. Multigrid methods overcome the critical slowing-down (CSD) phenomenon occuring with basic real-space relaxation methods. Several approaches employing the multigrid idea have appeared during recent years.

From the different multigrid methods available for the solution of the Schrödinger equation, we have picked up the Rayleigh Quotient Multigrid (RQMG) method introduced by Mandel and McCormick . This approach differs from full-approximation-storage (FAS) methods, as well as from those methods, where the eigenproblem is linearized.

In the RQMG method the coarse grid relaxation passes are performed so that the Rayleigh quotient calculated on the fine grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid and the coarse grid representation problem is avoided. Mandel and McCormick introduced the method for the solution of the eigenpair corresponding to the lowest eigenvalue. We have generalized it to the simultaneous solution of a desired number of lowest eigenenergy states by developing a scheme which keeps the eigenstates separated by the use of a penalty functional.

## II Numerical Methods

In our RQMG application the coarse grid relaxations are performed by the so-called coordinate relaxation method. It solves the discretized eigenproblem

$$Hu=\lambda Bu$$
(1)
by minimizing the Rayleigh quotient

$$\frac{u|H|u}{u|B|u}.$$
(2)
Above, $`H`$ and $`B`$ are matrix operators chosen so that the Schrödinger equation discretized on a real-space point grid with spacing $`h`$ is satisfied to a chosen order $`O(h^n)`$. In Eq. (2) $`u`$ is a vector containing the wave function values at the grid points. In the relaxation method, the current estimate $`u`$ is replaced by $`u^{}=u+\alpha d`$, where the search vector $`d`$ is simply chosen to be unity in one grid point and to vanish in all other points, and $`\alpha `$ is chosen to minimize the Rayleigh quotient. This leads to a simple ***For the sake of simplicity, the wave-functions, and thus $`\alpha `$, are here assumed real. We have implemented the complex case as well. quadratic equation for $`\alpha `$. A complete coordinate relaxation pass is then obtained by performing the minimization at each point in turn and these passes can be repeated until the lowest state is found with desired accuracy.

Naturally, also the coordinate relaxation suffers from CSD because of the use of local information only in updating $`u`$ in a certain point. In order to avoid it one applies the multigrid idea. In the multigrid scheme by Mandel and McCormick the crucial point is that coarse grid coordinate relaxation passes are performed so that the Rayleigh quotient calculated on the fine grid will be minimized. In this way there is no requirement for the solution to be well represented on a coarse grid. In practice, a coarse grid search substitutes the fine grid solution by

$$u_f^{}=u_f+\alpha I_c^fe_c,$$
(3)
where the subscripts $`f`$ and $`c`$ stand for the fine and coarse grids, respectively, and $`I_c^f`$ a prolongation operator interpolating the coarse grid vector to the fine grid. The Rayleigh quotient to be minimized is then

$`{\displaystyle \frac{u_f+\alpha I_c^fd_c|H_f|u_f+\alpha I_c^fd_c}{u_f+\alpha I_c^fd_c|B_f|u_f+\alpha I_c^fd_c}}=`$ (4)
$`{\displaystyle \frac{u_f|H_fu_f+2\alpha I_f^cH_fu_f|d_c+\alpha ^2d_c|H_cd_c}{u_f|B_fu_f+2\alpha I_f^cB_fu_f|d_c+\alpha ^2d_c|B_cd_c}}.`$ (5)
The second form is obtained by relating the coarse grid operators, $`H_c`$ and $`B_c`$, with the fine grid ones, $`H_f`$ and $`B_f`$, by the Galerkin condition

$$H_c=I_f^cH_fI_c^f;B_c=I_f^cB_fI_c^f;I_f^c=\left(I_c^f\right)^T.$$
(6)
The key point to note is that when $`H_fu_f`$ and $`B_fu_f`$ are provided from the fine grid to the coarse grid, the remaining integrals can be calculated on the coarse grid itself. Thus one really applies coordinate relaxation on the coarse grids to minimize the fine level Rayleigh quotient. This is a major departure from the earlier methods, which to some extent rely on the ability to represent the solution of some coarse grid equation on the coarse grid itself. Here, on the other hand, one can calculate the exact change in the Rayleigh quotient due to any coarse grid change, no matter how coarse the grid itself is. There is no equation whose solution would have to be representable.

In the MIKA package we have generalized the RQMG method to the simultaneous solution of several mutually orthogonal eigenpairs. The separation of the different states is divided into two or three subtasks. First, in order to make the coarse grid relaxations converge towards the desired state we apply a penalty functional scheme. Given the current approximations for the $`k`$ lowest eigenfunctions, the next lowest, $`(k+1)`$’th state is updated by minimizing the functional

$$\frac{u_{k+1}|H|u_{k+1}}{u_{k+1}|B|u_{k+1}}+\underset{i=1}{\overset{k}{}}q_i\frac{\left|u_i|u_{k+1}\right|^2}{u_i|u_iu_{k+1}|u_{k+1}}.$$
(7)
The minimization of this functional is equivalent to imposing the orthonormality constraints against the lower $`k`$ states, when $`q_i\mathrm{}`$. By increasing the shifts $`q_i`$ any desired accuracy can be obtained, but in order to obtain a computationally efficient algorithm a reasonable finite value should be used, for example

$$q_i=(\lambda _{k+1}\lambda _i)+\mathrm{Q},$$
(8)
where $`Q`$ is a sufficiently large positive constant. In our test calculations $`Q`$ is of the order of $`Q=0.5\mathrm{}2`$ Ha.

The substitution (3) is introduced in the functional (7) and the minimization with respect to $`\alpha `$ leads again to a quadratic equation. This time the coefficients contain terms due to the penalty part.

While the penalty functional keeps the states separated on the coarse levels, we apply a simple relaxation method (Gauss-Seidel) on the finest level. The Gauss-Seidel method converges to the nearest eigenvalue, so ideally no additional orthogonalizations would be needed. In practice, however, we use Gramm-Schmidt orthogonalizations and subspace rotations. However, the number of fine grid orthogonalizations remains quite plausible, for example, in comparison with the conjugate gradient search of eigenpairs employing only the finest grid .

The Kohn-Sham equations have to be solved self-consistently, i.e. the wave functions solved from the single-particle equation determine via the density (solution of the Poisson equation and the calculation of the exchange-correlation potential) the effective potential for which they should again be solved. To approach this self-consistency requires an optimized strategy so that numerical accuracy of the wave functions and the potential increase in balance, enabling the most efficient convergence . Our strategy in MIKA for self-consistency iterations is illustrated in Fig. 1. The Poisson equation for the Coulomb potential is solved also by the multigrid method.

## III Examples

We have demonstrated the performance of the MIKA scheme in calculating the electronic structures of small molecules and solid-state systems described by pseudopotentials. As a typical application Fig. 2 shows the electron density of the so-called deep state localized at a neutral, ideal vacancy in bulk Si. It was shown, that the accuracy of 1 meV for the total energy was reached after three or four V-cycles, and that the amount of cpu-time needed was of the same order as when applying state-of-the-art plane-wave codes. We obtained an average convergence rate of approximately one decade per self-consistency iteration. This is of the same order as those reported by Wang and Beck  in their FAS scheme or by Kresse and Furthmüller  in their plane-wave scheme employing self-consistency iterations. The convergence rate of one decade per self-consistency iteration is better than that obtained by Ancilotto et al.  in the FMG scheme and much better than the rate reached in the linearized multigrid scheme by Briggs et al. .

We have applied the RQMG method also for the calculation of positron states in solids. Fig. 3 shows how the delocalized positron state in the perfect $`\alpha `$-quartz is trapped in to a Si-vacancy. Positron states are a particularly simple case for our method, because only the lowest energy wave function needs to be calculated in a given potential, so that no orthogonalizations or penalty functionals are needed. Moreover, in a simple scheme an electron density calculated without the influence of the positron can be used as the starting point . However, even for the positron states the superior performance of the multigrid method in comparison with straightforward relaxation schemes is evident. For example, we have calculated the positron state in the Si-vacancy in bulk Si using a supercell containing 1727 atoms. The solution of the positron wave-function using the RQMG-method took less than a minute of cpu-time on a typical work station. To put this in proper context, J. E. Pask et al  report a similar calculation, based on the finite element method but without multigrid acceleration, for a supercell containing 4096 Cu atoms. The result converged within 1 ps took ’ just 14.3 hr ’ of CPU time.

We have also applied the MIKA scheme in two-dimensional problems for quantum dots employing the current-spin-density functional theory (CSDFT), see Ref. . Moreover, we have implemented the RQMG-method in cylindrical coordinates enabling very efficient and accurate calculations for atomic chains, or systems which can be described using axisymmetric jellium models. Fig. 4 shows a selected wavefunction of a system where a chain of four carbon atoms is sandwhiched between two planar jellium leads.

## IV Summary and outlook

In the MIKA program package the RQMG method introduced by Mandel and McCormick  is generalized for the simultaneous solution of a desired number of lowest eigenenergy states. The approach can be viewed to belong to a third group of multigrid methods, in addition to FAS and techniques where the eigenproblem is linearized. In principle, one can use arbitrarily coarse grids in RQMG, whereas in the other multigrid methods one has to be able to represent all the states also on the coarsest grid.

We are convinced that our method will compete with the standard plane-wave methods for electronic structure calculations. However, some straightforward programming is still required. Implementation of the Hellmann-Feynman forces, required for the optimization of the ionic structures is under way.

During the RQMG V-cycle, the states are all relaxed simultaneously and independently of each other. A parallelization over states would therefore be natural to implement on a shared memory architecture. We have parallelized the MIKA codes over k-points, and over real-space domains. The domain decomposition is the appropriate method for distributed memory parallel computers.

###### Acknowledgements.

We acknowledge the contributions by Henri Saarikoski, Paula Havu, Esa Räsänen, Tero Hakala, and Sampsa Riikonen in sharing their experience of the use of the MIKA package in different applications and preparing the figures 3 (T.H.) and 4 (P.H.). T.T. acknowledges financial support by the Vilho, Yrjö and Kalle Väisälä foxundation. This research has been supported by the Academy of Finland through its Centre of Excellence Programme (2000 - 2005).


text:
# References

hep-ph/0205032

A Full Determination of the Neutrino Mass Spectrum
from Two-zero Textures of the Neutrino Mass Matrix

Zhi-zhong Xing ***Electronic address: xingzz@mail.ihep.ac.cn

CCAST (World Laboratory), P.O. Box 8730, Beijing 100080, China

and Institute of High Energy Physics, Chinese Academy of Sciences,
P.O. Box 918 (4), Beijing 100039, China Mailing address

## Abstract

We show that it is possible to fully determine the neutrino mass spectrum from two-zero textures of the neutrino mass matrix. As a consequence, definite predictions can be obtained for the neutrinoless double beta decay.

In the flavor basis where the charged lepton mass matrix is diagonal, a phenomenological analysis of two-zero textures and Majorana CP-violating phases of the neutrino mass matrix has recently been carried out by the author . Seven of such textures were found to be in accord with current experimental data on atmospheric , solar  and reactor  neutrino oscillations. A classification of two-zero textures of the neutrino mass matrix was originally done by Frampton, Glashow and Marfatia . In two follow-up works, Frampton, Oh and Yoshioka  have proposed a model with three Higgs triplets to gain texture zeros of the neutrino mass matrix; and Kageyama, Kaneko, Shimoyama and Tanimoto  have incorporated those phenomenologically acceptable textures with the seesaw mechanism.

As an addendum to Ref. , the present letter aims to show that it is actually possible to fully determine the neutrino mass spectrum from two-zero textures of the neutrino mass matrix. This important point was not observed in Ref.  and Refs. , where only the ratios of three neutrino masses were calculated. Given the neutrino mass spectrum, a definite prediction for the neutrinoless double beta decay can be made. Taking account of current experimental data, which remain rather rough, we illustrate the typical magnitudes of neutrino masses and that of the neutrinoless double beta decay for each of the seven two-zero textures of the neutrino mass matrix.

The Majorana neutrino mass matrix can be written as

$$M=V\left(\begin{array}{ccc}m_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right)V^\mathrm{T}$$
(1)
in the flavor basis where the charged lepton mass matrix is diagonal, where $`m_i`$ (for $`i=1,2,3`$) are physical neutrino masses and $`V`$ stands for the lepton flavor mixing matrix linking the neutrino mass eigenstates $`(\nu _1,\nu _2,\nu _3)`$ to the neutrino flavor eigenstates $`(\nu _e,\nu _\mu ,\nu _\tau )`$. Following Ref. , we parametrize $`V`$ as

$$V=\left(\begin{array}{ccc}c_xc_z& s_xc_z& s_z\\ c_xs_ys_zs_xc_ye^{i\delta }& s_xs_ys_z+c_xc_ye^{i\delta }& s_yc_z\\ c_xc_ys_z+s_xs_ye^{i\delta }& s_xc_ys_zc_xs_ye^{i\delta }& c_yc_z\end{array}\right)\left(\begin{array}{ccc}e^{i\rho }& 0& 0\\ 0& e^{i\sigma }& 0\\ 0& 0& 1\end{array}\right),$$
(2)
where $`s_x\mathrm{sin}\theta _x`$, $`c_y\mathrm{cos}\theta _y`$, and so on. In this “standard” parametrization of $`V`$, the Dirac phase $`\delta `$ controls the strength of CP violation in normal neutrino oscillations, while the Majorana phases $`\rho `$ and $`\sigma `$ are relevant to the neutrinoless double beta decay . Without loss of generality, three mixing angles ($`\theta _x,\theta _y,\theta _z`$) can all be arranged to lie in the first quadrant. Arbitrary values are possible for three CP-violating phases ($`\delta ,\rho ,\sigma `$). Note that $`M`$ is symmetric, thus it consists totally of six independent complex entries. Assuming two of them to vanish, one may figure out fifteen distinct textures of $`M`$. Only seven textures, as listed in Table 1, were found to be compatible with current experimental data and empirical hypotheses . Their predictions for the ratios of three neutrino masses,

$$\xi \frac{m_1}{m_3},\zeta \frac{m_2}{m_3},$$
(3)
are quoted from Ref.  and listed in the same table. Our present purpose is to determine $`m_1`$, $`m_2`$ and $`m_3`$ separately, so as to fully fix the neutrino mass spectrum.

Given $`\xi `$ and $`\zeta `$, the absolute values of three neutrino masses can be extracted from the measured mass-squared differences of solar and atmospheric neutrino oscillations (i.e., $`\mathrm{\Delta }m_{\mathrm{sun}}^2`$ and $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$). Because current Super-Kamiokande and SNO data strongly support the hypothesis that solar and atmospheric neutrino oscillations are dominated respectively by $`\nu _e\nu _\mu `$ and $`\nu _\mu \nu _\tau `$ transitions , one may simply define $`\mathrm{\Delta }m_{\mathrm{sun}}^2`$ and $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ as

$`\mathrm{\Delta }m_{\mathrm{sun}}^2`$ $``$ $`\left|m_2^2m_1^2\right|,`$
$`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ $``$ $`\left|m_3^2m_2^2\right|.`$ (4)
In Refs. , the large-angle MSW solution to the solar neutrino problem has been taken into account. Thus $`R_\nu \mathrm{\Delta }m_{\mathrm{sun}}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2𝒪(10^2)`$ has been used as a crucial criterion to single out the textures of $`M`$ which are phenomenologically favored. From Eqs. (3) and (4), it is straightforward to obtain

$$m_3=\sqrt{\frac{\mathrm{\Delta }m_{\mathrm{atm}}^2}{|\zeta ^21|}},$$
$`(5\mathrm{a})`$
or equivalently

$$m_3=\sqrt{\frac{\mathrm{\Delta }m_{\mathrm{sun}}^2}{|\xi ^2\zeta ^2|}}.$$
$`(5\mathrm{b})`$
Once $`m_3`$ is determined from Eq. (5a) or (5b), the magnitudes of $`m_1`$ and $`m_2`$ can be fixed by use of Eq. (3) or with the help of the formulas presented in Ref. .

To be specific, let us calculate $`m_3`$ for each of the seven textures of $`M`$ listed in Table 1. We use $`\mathrm{\Delta }m_{\mathrm{atm}}^23\times 10^3\mathrm{eV}^2`$  and $`\mathrm{\Delta }m_{\mathrm{sun}}^25\times 10^5\mathrm{eV}^2`$  as typical inputs. Note that the number of $`\mathrm{\Delta }m_{\mathrm{sun}}^2`$ results from the best global fit of present solar neutrino oscillation data (including the latest SNO neutral current measurement ) in the large-angle MSW mechanism, and the corresponding value of $`\theta _x`$ is $`\theta _x30^{}`$. Furthermore, $`\theta _y40^{}`$ and $`\theta _z5^{}`$ are typically taken, although the latter may be much smaller. We assume the unknown CP-violating phase $`\delta `$ to be around $`90^{}`$ in most cases. As the uncertainty associated with $`\mathrm{\Delta }m_{\mathrm{atm}}^2`$ is expected to be smaller than that associated with $`\mathrm{\Delta }m_{\mathrm{sun}}^2`$, we choose to use Eq. (5a) instead of Eq. (5b) in the calculation of $`m_3`$.

Pattern $`\mathrm{A}_1:`$ Because of $`s_z1`$ and $`t_xt_y`$, $`\zeta ^21`$ is naturally anticipated. Therefore we obtain

$$m_3\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}\mathrm{\hspace{0.33em}5.5}\times 10^2\mathrm{eV}.$$
(6)
Using the typical inputs $`\theta _x=30^{}`$, $`\theta _y=40^{}`$ and $`\theta _z=5^{}`$, we arrive at

$$m_1\mathrm{\hspace{0.33em}2.3}\times 10^3\mathrm{eV},m_2\mathrm{\hspace{0.33em}7.0}\times 10^3\mathrm{eV}.$$
(7)
This quasi-hierarchical spectrum of neutrino masses implies that it is in practice impossible to detect the neutrinoless double beta decay.

Pattern $`\mathrm{A}_2`$: Once again $`\zeta ^21`$ holds. Then we obtain the same value of $`m_3`$ as that given in Eq. (6). Using the same inputs as above, we find

$$m_1\mathrm{\hspace{0.33em}3.3}\times 10^3\mathrm{eV},m_2\mathrm{\hspace{0.33em}1.0}\times 10^2\mathrm{eV}.$$
(8)
One can see that the neutrino mass spectra of Patterns $`\mathrm{A}_1`$ and $`\mathrm{A}_2`$ are quite similar.

Pattern $`\mathrm{B}_1`$: With the help of Eq. (5a), we obtain

$$m_3\sqrt{\frac{\mathrm{\Delta }m_{\mathrm{atm}}^2}{|t_y^41|}}\mathrm{\hspace{0.33em}7.7}\times 10^2\mathrm{eV},$$
(9)
where $`\theta _y=40^{}`$ has been used. In the lowest-order approximation, we have

$$m_1m_2\mathrm{\hspace{0.33em}5.4}\times 10^2\mathrm{eV}.$$
(10)
This quasi-degenerate spectrum of neutrino masses might have useful hints at possible flavor symmetries and their breaking schemes, which are expected to be responsible for the generation of lepton masses  and associated with the origin of leptogenesis . It is also worth mentioning that the values of $`m_i`$ in Eqs. (9) and (10) are compatible with the present direct-mass-search experiments , in particular for the electron neutrino.

Pattern $`\mathrm{B}_2`$: In analogy to Eq. (9), $`m_3`$ reads as

$$m_3t_y^2\sqrt{\frac{\mathrm{\Delta }m_{\mathrm{atm}}^2}{|t_y^41|}}\mathrm{\hspace{0.33em}5.4}\times 10^2\mathrm{eV},$$
(11)
where $`\theta _y=40^{}`$ has been used. To lowest order, we obtain

$$m_1m_2\mathrm{\hspace{0.33em}7.7}\times 10^2\mathrm{eV}.$$
(12)
We see that the neutrino mass spectra of Patterns $`\mathrm{B}_1`$ and $`\mathrm{B}_2`$ have much similarity too.

Pattern $`\mathrm{B}_3`$: To lowest order, the neutrino mass spectrum in this pattern is identical to that in Pattern $`\mathrm{B}_1`$.

Pattern $`\mathrm{B}_4`$: To lowest order, the neutrino mass spectrum in this pattern is identical to that in Pattern $`\mathrm{B}_2`$. We have seen that the phenomenological consequences of Patterns $`\mathrm{B}_1`$, $`\mathrm{B}_2`$, $`\mathrm{B}_3`$ and $`\mathrm{B}_4`$ are nearly the same . It is particularly difficult to distinguish between Patterns $`\mathrm{B}_1`$ and $`\mathrm{B}_3`$, or between Patterns $`\mathrm{B}_2`$ and $`\mathrm{B}_4`$.

Pattern C: In the lowest-order approximation, we obtain

$$m_3|t_{2y}|s_z\sqrt{\frac{\mathrm{\Delta }m_{\mathrm{atm}}^2}{t_x|t_x+2t_{2y}s_zc_\delta |}}.$$
(13)
One can see that the magnitude of $`m_3`$ depends not only upon three mixing angles $`(\theta _x,\theta _y,\theta _z)`$ but also upon the Dirac phase $`\delta `$. Therefore a careful analysis of the parameter space is required for Pattern C , so as to fit current experimental data and to fix the allowed range of $`m_3`$. The values of $`\delta `$ and $`\theta _z`$ are unfortunately unknown at present. Just for the purpose of illustration, we typically take $`\theta _x=\theta _y=44.8^{}`$, $`\theta _z=5^{}`$ and $`\delta =90^{}`$. In this special case, one may get $`R_\nu 0.03`$ – a correct order of $`\mathrm{\Delta }m_{\mathrm{sun}}^2/\mathrm{\Delta }m_{\mathrm{atm}}^2`$ . Then we arrive at

$$m_3\frac{t_{2y}s_z}{t_x}\sqrt{\mathrm{\Delta }m_{\mathrm{atm}}^2}\mathrm{\hspace{0.33em}0.7}\mathrm{eV}.$$
(14)
Using the same inputs, we get $`m_1m_2m_30.7\mathrm{eV}`$ to a high degree of accuracy. This result indicates that three neutrino masses are essentially degenerate, and their magnitude can be of $`𝒪(1)`$ eV. Therefore it is rather sensitive to the neutrinoless double beta decay.

Indeed two-zero textures of the neutrino mass matrix allow us to obtain definite predictions for the neutrinoless double beta decay, whose effective mass term is a simple function of neutrino masses and flavor mixing parameters:

$$m_{ee}=m_3\left|V_{e1}^2\xi +V_{e2}^2\zeta +V_{e3}^2\right|.$$
(15)
Using the parametrization of $`V`$ in Eq. (2) and the expressions of $`\xi `$ and $`\zeta `$ in Table 1, we can calculate $`m_{ee}`$ for each of the seven patterns of $`M`$. The approximate analytical results are listed in Table 2. Some comments are in order.

(1) $`m_{ee}0`$ holds for Patterns $`\mathrm{A}_1`$ and $`\mathrm{A}_2`$. This is obviously true, as $`M_{ee}=0`$ has been taken in both patterns.

(2) The sizes of $`m_{ee}`$ in Patterns $`\mathrm{B}_1`$ and $`\mathrm{B}_3`$ are essentially identical: $`m_{ee}m_1m_25.4\times 10^2`$ eV for $`\theta _y=40^{}`$. So are the sizes of $`m_{ee}`$ in Patterns $`\mathrm{B}_2`$ and $`\mathrm{B}_4`$: $`m_{ee}m_1m_27.7\times 10^2`$ eV for $`\theta _y=40^{}`$.

(3) If $`\theta _x=\theta _y=44.8^{}`$, $`\theta _z=5^{}`$ and $`\delta =90^{}`$ are typically taken, one will arrive at $`m_{ee}m_1m_2m_30.7`$ eV for Pattern C At this point we notice that there is a typing error associated with $`|M_{ee}|`$ in Eq. (29) of Ref. . The correct result should be $`|M_{ee}|m_3\sqrt{14c_\delta /(t_{2x}t_{2y}s_z)+4/(t_{2x}^2t_{2y}^2s_z^2)}`$.. This result is apparently consistent with the alleged evidence for the neutrinoless double beta decay , $`0.05`$ eV $`m_{ee}`$ 0.84 eV, at the $`95\%`$ confidence level. Nevertheless, it might be premature to take this measurement too seriously . To be more conservative, we take account of the relatively reliable experimental upper bound $`m_{ee}<0.35`$ eV (at the $`90\%`$ confidence level ) and find that the typical result $`m_{ee}=0.7`$ eV is actually disfavored. It should be noted, however, that there does exist the appropriate parameter space for Pattern C , in which both $`R_\nu 𝒪(10^2)`$ and $`m_{ee}<0.35`$ eV can be satisfied.

As pointed out in Ref. , the seven patterns of $`M`$ can be classified into three distinct categories: A (with $`\mathrm{A}_1`$ and $`\mathrm{A}_2`$), B (with $`\mathrm{B}_1`$, $`\mathrm{B}_2`$, $`\mathrm{B}_3`$ and $`\mathrm{B}_4`$), and C. A definite measurement of the neutrinoless double beta decay at the level $`m_{ee}𝒪(0.01)`$ eV to $`𝒪(0.1)`$ eV will rule out category A. On the other hand, more precise data of neutrino oscillations will help reduce the uncertainties associated with $`m_{ee}`$ predicted in categories B and C. It might even be possible to distinguish between categories B and C (and therefore discard one of them), if the magnitude of $`m_{ee}`$ is experimentally determined in the future.

In summary, we have shown that a full determination of the neutrino mass spectrum is indeed possible from two-zero textures of the neutrino mass matrix $`M`$. This important observation indicates that two-zero textures of $`M`$ have much more predictability than previously expected. In particular, one can get definite predictions for the neutrinoless double beta decay. We hope that a variety of proposed precision measurements of neutrino oscillations and lepton-number-violating processes may finally allow us to pin down the unique texture of lepton mass matrices.

The author would like to thank H.J. He, D. Marfatia and T. Yoshikawa for useful discussions. He is also grateful to H.V. Klapdor-Kleingrothaus for pointing out an error in the list of references. This work was supported in part by the National Natural Science Foundation of China.


text:
# Observation of minority spin character of the new electron doped manganite La0.7Ce0.3MnO3 from tunneling magnetoresistance

## Abstract

We report the magnetotransport characteristics of a trilayer ferromagnetic tunnel junction build of an electron doped manganite (La0.7Ce0.3MnO3) and a hole doped manganite (La0.7Ca0.3MnO3). At low temperatures the junction exhibits a large positive tunneling magnetoresistance (TMR), irrespective of the bias voltage. At intermediate temperatures below $`T_C`$ the sign of the TMR is dependent on the bias voltage across the junction. The magnetoresistive characteristics of the junction strongly suggest that La0.7Ce0.3MnO3 is a minority spin carrier ferromagnet with a high degree of spin polarization, i.e. a transport half metal.

There has been a lot of interest recently in the hole doped rare-earth manganites, where the rare-earth in the insulating parent compound is partially replaced by a divalent cation (such as Ca, Ba, Sr, Pb etc.) Coey et al. (1999). Around 30% hole doping most of these compounds such as La0.7Ca0.3MnO3, La0.7Ba0.3MnO3 have a half-metallic ferromagnetic ground state with a large magnetoresistance (MR) associated with the ferromagnetic transition temperature. Conversely, doping electrons by substituting the rare-earth atom by tetravalent Ce also drives the system into a ferromagnetic metallic ground state Mandal and Das (1997); Raychaudhuri et al. (1999); Mitra et al. (2001a). Unlike the hole doped compound where manganese is in a mixture of Mn3+ and Mn4+ valence states the tetravalent Ce drives the compound in a mixture of Mn3+ and Mn2+ valencies. The presence of Mn2+/Mn3+ valencies induced by electron doping in La0.7Ce0.3MnO3 was confirmed recently Mitra et al. (a). There exists an inherent symmetry between Mn4+ and Mn2+ as both are non Jahn-Teller ions, whereas Mn3+ is a Jahn-Teller ion. Moreover, La0.7Ca0.3MnO3 and La0.7Ce0.3MnO3 both have a Curie temperature of $`T_C`$ 250 K.

Beyond the phenomenon of colossal MR, novel properties arising from the interplay of spin, charge and orbital coupling and the competition of closely related energy scales make the manganites potential candidates for novel electronic devices. Devices exhibiting both large positive and negative MR have been fabricated by integrating doped manganites with other oxide ferromagnets. More recently rectifying characteristics could be demonstrated in La0.7Ce0.3MnO3/SrTiO3(STO)/La0.7Ca0.3MnO3 tunnel junctions Mitra et al. (2001b) and Nb-STO/(La,Ba)MnO3 junctions Tanaka et al. (2002). However, a detailed understanding of the spin dependent electronic structure of these materials is important to exploit the complete potential of these compounds.

In this paper we report the magnetotransport properties of a La0.7Ce0.3MnO3/STO($`50`$Å)/La0.7Ca0.3MnO3 tunnel junction grown onto conducting 0.5% Nb doped SrTiO3 (Nb-STO) single crystalline substrate through pulsed laser deposition, whose fabrication has been reported earlier Mitra et al. (2001b). The device structure is illustrated in the inset of Fig. 1(a). The Nb-STO substrate acts as a conducting underlayer allowing us to measure the transport properties of the tunnel junction with current perpendicular to plane (CPP) geometry. In contrast to the hole doped compounds not much work Min et al. (2001); Kang et al. (2001) could so far be done to understand the electronic structure of La0.7Ce0.3MnO3. A major experimental obstacle stems from the fact that La0.7Ce0.3MnO3 has so far been synthesized in single phase only in thin film form through high-energy pulsed laser deposition Mitra et al. (2001a). This is a restriction on many of the conventional probes and calls for alternative ways to extract information regarding the spin polarized electronic structure in these materials. Towards this end we have used magnetotransport properties of tunnel junctions comprising of two ferromagnetic electrodes separated by a thin insulating tunnel barrier. Most importantly, we wanted to establish whether this electron doped ferromagnet has minority spin carriers (MISC) (where the magnetization is antiparallel to the spin) or majority spin carriers (MASC) (magnetization parallel to the spin) at the Fermi energy $`E_F`$. The resistance across such a ferromagnetic tunnel junction (FTJ) depends on the relative orientation of the spins at $`E_F`$ in the two electrodes. In zero field the magnetizations of the two ferromagnetic electrodes in an ideal FTJ are predominantly antiparallel owing to magnetostatic interaction between the ferromagnetic layers with in-plane magnetization. With the application of a magnetic field the magnetizations become parallel. Thus, a tunnel junction where both electrodes are either MISC or MASC ferromagnets exhibits negative tunneling magnetoresistance (TMR). This is observed in manganite tunnel junctions of identical materials separated by a thin STO insulating barrier Gupta et al. (1996); Sun et al. (1996). In contrast, for a tunnel junction with one MISC and one MASC ferromagnetic electrode the MR is positive. This has experimentally been observed in La0.7Sr0.3MnO3/STO/Fe3O4 tunnel junctions Ghosh et al. (1998) showing large positive MR, where Fe3O4 is a MISC half-metallic ferrimagnet and La0.7Sr0.3MnO3 is the MASC half-metallic ferromagnet. Here, we investigate the spin character of La0.7Ce0.3MnO3 at $`E_F`$ from the magnetotransport properties of the tunnel junction using the MASC ferromagnet La0.7Ca0.3MnO3 as a spin analyzer.

Fig. 1(a) shows the magnetic field dependence of the current versus voltage ($`IV`$) curve across the tunnel junction measured in the CPP geometry at 300 K. We do not see a significant TMR, as expected, since at this temperature both La0.7Ce0.3MnO3 and La0.7Ca0.3MnO3 are paramagnetic semiconductors and the device behaves like a rectifying diode. However, in a field of 7.5 T the spin disorder scattering of a single layer is reduced, and the in-plane MR is quite large Coey et al. (1999). Fig. 1(b) shows the magnetic field dependence of the tunneling $`IV`$ curve taken at 100 K, in zero field and in a field of 2 T. The magnetic field dependence of the $`IV`$ curves (both positive and negative bias) clearly shows a bias dependent MR, i.e., below a threshold bias voltage (or current), we find a negative MR and above this a positive MR prevails. The field dependence of the $`IV`$ curves taken at 48 K \[Fig. 1(c)\] does not show any bias dependence, exhibiting a positive MR at all voltages. In order to rule out any significant influence of the Nb-STO/La0.7Ce0.3MnO3 bottom junction on the results of the tunnel junction we also deposited single layers of La0.7Ce0.3MnO3 (and La0.7Ca0.3MnO3) of identical thickness upon Nb-STO. As a main result, the MR of the Nb-STO/La0.7Ce0.3MnO3 bottom junction is several orders of magnitude smaller compared to that of the FTJ and is always negative. Moreover, the $`IV`$ curve of the former is highly asymmetric \[inset Fig. 1(c)\] down to the lowest temperature.

First we concentrate on the $`IV`$ characteristics of the tunnel junction at the lowest temperature (48 K). Here, the positive TMR is intriguing. In hole doped manganites, three of the Mn-3d electrons form the localized $`t_{2g}`$ band. The remaining electrons occupy the conducting $`e_g`$ band which is energetically higher than the $`t_{2g}`$ in a (nearly) cubic crystal field. The crystal field splitting energy is estimated to be $`\mathrm{\Delta }_{cf}`$ $``$ 1.8 eV Coey et al. (1999). The $`e_g`$ band splits further due to Jahn-Teller distortion into two sub-bands, $`e_g^1`$ and $`e_g^2`$, which are separated by the Jahn-Teller splitting energy, $`\delta _{JT}`$ $``$ 1.2 eV Coey et al. (1999); Satpathy et al. (1996). Jahn-Teller distortion also causes a splitting of the $`t_{2g}`$ band by $`\delta _{JT}^{}`$. Moreover, Hund’s rule coupling removes the spin degeneracy in the ferromagnetic state. The resulting separation of the spin-up ($`e_g`$) and spin-down ($`e_g`$) bands is denoted by $`U_H`$ whereas the separation of the $`t_{2g}`$ and $`t_{2g}`$ bands is labeled $`U_H^{}`$. The Mn-3d spin dependent density of states (DOS) of La0.7Ca0.3MnO3 at low temperature is sketched in Fig. 2(a). The bandwidth of the conducting $`e_g`$ sub-bands is of the order of 1 eV Coey et al. (1999). $`U_H`$ is estimated from band structure calculations to be of the order of 2 eV in the undoped compound—a value that may decrease with doping in the hole doped compound Okimoto et al. (1997). In hole doped manganites

like La0.7Ca0.3MnO3, which have 0.7 electrons in the $`e_g`$ band, the conduction electrons predominantly occupy the lowest sub-band $`e_g^1`$ and the hole doped manganite is consequently a MASC ferromagnet \[Fig. 2(a)\].

In La0.7Ce0.3MnO3, there is clear evidence for electron doping on the Mn-site in single layer epitaxial films, e.g., from x-ray absorption spectroscopy Mitra et al. (a). Hence, the $`e_g^1`$ sub-band is completely filled. For the remaining additional (doped) electrons two scenarios are possible: i) weak Hund’s rule coupling $`U_H^{}<\mathrm{\Delta }_{cf}+\delta _{JT}`$ and ii) strong Hund’s rule coupling $`U_H^{}>\mathrm{\Delta }_{cf}+\delta _{JT}`$. In the first case, $`t_{2g}`$ is energetically lower than $`e_g^2`$ and the former will get partially filled \[Fig. 2(b)\]. The compound will be a MISC ferromagnet in an intermediate spin state. In the second case, the remaining electrons will occupy the $`e_g^2`$ sub-band resulting in a MASC ferromagnet in high spin state \[Fig. 2(c)\]. The observation of a positive TMR in La0.7Ce0.3MnO3/STO/La0.7Ca0.3MnO3 at low temperature definitely favors the first scenario with antiparallel spins at $`E_F`$ for fields high enough to align the magnetizations within the two ferromagnetic layers \[cf. Fig. 2(a) and (b)\]. This result is supported by the known energy values given above if $`U_H^{}U_H`$ is assumed. Band structure calculations Satpathy et al. (1996) also predict the near-degeneracy of the energy positions of the $`e_g^2`$ and the $`t_{2g}`$ sub-bands.

The minority spin character and the related intermediate spin

state observed in La0.7Ce0.3MnO3 are intriguing for several reasons. First, due to the usually large on-site Hund’s rule coupling this is rarely observed in manganese compounds where Mn is in the divalent state including compounds such as MnO. Secondly, the transport and magnetic properties of the hole doped manganites are conventionally understood from this large on-site Hund’s rule coupling in manganese and the electron-lattice coupling arising from the Jahn-Teller effect in Mn3+. The ferromagnetic ground state in hole doped compounds like La0.7Ca0.3MnO3 is driven by the double exchange interaction between Mn3+ and Mn4+. This, however, is unlikely to be the dominant mechanism in La0.7Ce0.3MnO3 if the manganese is in intermediate spin state. The striking similarity between the magnetic and transport properties of La0.7Ce0.3MnO3 and its hole doped counterpart Raychaudhuri et al. (1999) raises questions anew regarding the origin of colossal MR in these compounds. It should be pointed out in this context that even in the hole doped compound the existence of Mn3+ and Mn4+ in the pure high spin state has been questioned recently. From point contact Andreev reflection measurements Nadgorny et al. (2001) the existence of minority spin carriers in the hole doped compound La0.7Sr0.3MnO3 was concluded. Studies on tunnel junctions using a superconducting electrode and La0.7Sr0.3MnO3 Worledge and Geballe (2000) also showed the polarization $`P`$ $`=(n_{}(E_F)n_{}(E_F))/(n_{}(E_F)+n_{}(E_F))`$ ($`n_{}`$ and $`n_{}`$ being the spin up and spin down DOS) to be of the order of 72%, which is much smaller than the polarisation expected for a half-metal Park et al. (1998). Existing band structure calculations Pickett and Singh (1996) predict an even lower spin polarization. These results indicate that the manganese is not in the pure high spin state in both electron and hole doped manganites with perovskite crystal structure. Therefore the general explanation of transport and magnetic properties of colossal MR manganites in terms of double exchange as the dominant mechanism needs a closer examination regardless of their high half-metallic transport character.

Fig. 3 shows the resistance versus field ($`RH`$) of the

tunnel junction at 100 K. At a very small bias current ($`I=0.06\mu `$A) the MR is negative while at a large current ($`I=3\mu `$A) the MR is positive. A description of the bias dependence of the MR at 100 K \[Fig. 1(b)\] and 200 K (not shown) is beyond the ground state picture described in Fig. 2. Here, a detailed knowledge of the band structure and the temperature evolution of the bands in La0.7Ce0.3MnO3 would be required. Since at high bias the large electric field across the insulating tunnel barrier shifts their Fermi energies considerably with respect to each other, we first concentrate on the temperature dependent TMR observed at low bias. The crossover of the TMR from positive to negative at low bias suggests a change in the dominant spin character of the ferromagnet at $`E_F`$ with temperature. This could happen for two reasons. First, in transition metal oxides the interplay between competing energy scales can induce a crossover from one spin state to another with temperature. For example, in LaCoO3 where the crystal field energy slightly exceeds the Hund’s rule coupling energy, Co is in a low spin, effectively non-magnetic state at low temperature. Around 100 K the cobalt ions undergo a crossover from a low spin state to an intermediate/high spin state accompanied by a change in magnetic and electric properties Imada et al. (1998). Since the $`t_{2g}`$ and $`e_g^2`$ bands are at similar energy levels (and may even overlap) in the manganites, a temperature dependent crossover from an intermediate spin state to a high spin state cannot be ruled out in La0.7Ce0.3MnO3. It may be caused, e.g., by decreasing $`\mathrm{\Delta }_{cf}`$ with increasing temperature. This needs to be substantiated through other measurements. Secondly, even in the absence of such a crossover in the Mn spin state the dominant spin character of the carriers at $`E_F`$ in a ferromagnet can change. For example, in a Stoner ferromagnet the exchange splitting between spin up and spin down bands decreases with increasing temperature, which changes the filling of the up and down spin bands. Therefore, subtle features in the DOS as a function of energy $`E`$ can modify the relative ratio of up ($`e_g`$) and down ($`t_{2g}`$) spins at $`E_F`$ without changing the overall spin state of the underlying Mn-ions.

We now focus on the bias dependence of the MR at 100 K and 200 K. A bias voltage $`V_b`$ on a metallic tunnel junction shifts the Fermi levels of the two electrodes by $`eV_b`$. The tunneling of electrons across an insulating barrier, however, occurs at equal energy levels. Hence, the tunneling probabilities for the magnetization of the two electrodes parallel (antiparallel) to each other, $`T_{}`$ ($`T_{}`$), are given by

$`T_{}`$ $``$ $`n_{}(E^+)n_{}^{}(E^{})+n_{}(E^+)n_{}^{}(E^{}),`$
$`T_{}`$ $``$ $`n_{}(E^+)n_{}^{}(E^{})+n_{}(E^+)n_{}^{}(E^{})`$
respectively, where $`E^+=E_F+\frac{eV_b}{2}`$ and $`E^{}=E_F\frac{eV_b}{2}`$. Here, $`n_{}(E)`$ and $`n_{}(E)`$ are the spin up and spin down DOS in La0.7Ca0.3MnO3 and $`n_{}^{}(E)`$ and $`n_{}^{}(E)`$ are the corresponding ones in La0.7Ce0.3MnO3. The relative magnitude of $`T_{}`$ and $`T_{}`$ depend on the details of the spin up and spin down DOS. At 100 K and 200 K, the negative MR at low bias voltage suggests that $`T_{}>T_{}`$ giving a negative MR. However, at the same temperatures, at large $`V_b`$ the MR may change sign if this inequality is altered owing to features in the energy dependent DOS. Moreover, small changes in the relative energies of the $`t_{2g}`$ and $`e_g^2`$ bands due to $`V_b`$ may change their occupancy which may result in a bias dependent crossover. Here, a detailed understanding can emerge only when the DOS in La0.7Ce0.3MnO3 is known in sufficient detail.

In summary, we have reported the magnetotransport properties of a tunnel junction made of electron and hole doped manganites. The observed minority spin carrier transport in La0.7Ce0.3MnO3 is of fundamental interest to understand the interplay of the Hund’s rule coupling energy with other energy scales such as Jahn-Teller energy and crystal field energy in doped manganites. We believe that this result will motivate further studies of the exact electronic structure of both electron and hole doped manganites and may open up an alternative approach towards spintronics.

The authors would like to acknowledge Andy Mackenzie, F. Steglich, S.K. Dhar and R. Pinto for helpful discussions and encouragement. P.R. and C.M. thank the Leverhulme Trust and DFG (Grant No. SFB 422), respectively, for financial support.


text:
# THE LOWER MAIN SEQUENCE AND MASS FUNCTION OF THE GLOBULAR CLUSTER MESSIER 41

## 1 Introduction

Analysis of the Hubble Deep Field (HDF) (Williams et al. 1995), which is to date the deepest optical image ever taken, has produced a remarkable series of important results both in Galactic and extragalactic astronomy (see Ferguson et al. 2000 for a summary). With the original HDF as a model, we proposed to carry out an HDF–like project in the nearest globular cluster, M4. We were successful in obtaining a roughly similar amount of telescope time (123 orbits) for the 2 filters F606W (V – 127,400 sec) and F814W (I – 192,400 sec) as had been granted for the HDF. Our plan was to reach $`\mathrm{V}30`$, $`\mathrm{I}28`$ with s/n about 3 in each filter in a single M4 field located 6 core radii ($`5\mathrm{}`$) from the cluster center – the same field which we observed in 1995 (GO 5461).

The main science driver was to search for the oldest and coolest white dwarfs (WDs) in M4. The discussion of the WD population in M4 is contained in the following Letter (Hansen et al. 2002). Other science goals are to attempt to identify the termination of the hydrogen burning MS and hence set the location of the brown dwarf boundary for metal–poor stars, establish the MF in the cluster down to very low masses and use the cluster WDs to extend the MF beyond the cluster turnoff, investigate the binary frequency in the cluster both through the location of stars in the CMD and via variability, explore for planets through stellar occultations, and examine the inner halo populations (Rich et al. 2002, in preparation) and compare it with that of the cluster.

## 2 The Data and Its Reduction

The complete data set consists of 15$`\times `$2600 sec exposures in F555W, 9$`\times `$800 sec in F814W (1995 – part of GO 5461 from cycle 4 – see Richer et al. 1995, 1997; Ibata et al. 1999) together with 98$`\times `$1300 sec in F606W and 148$`\times `$1300 sec obtained in F814W (2001 – GO 8679 in cycle 9). All images were secured with the same roll angle in the same field of M4 located at about 6 core radii from the cluster center. We also processed and reduced the small set of F814W images from GO 8153 but in the final reductions they were not used.

The new images as well as the earlier frames were preprocessed according to the recipes given in Stetson (1987) and Stetson et al. 1998. The individual frames were transformed to the coordinate system of the first new epoch frame using the MONTAGE2 program. Once the frames were all on the same coordinate system, those corresponding to a particular epoch and filter were averaged together with pixel rejection – the n highest pixels being rejected to eliminate cosmic ray contamination. From the HST WFPC2 Manual the mean number of pixels on a given chip affected by a cosmic ray hit in an 1800 sec exposure is about 20,000. Hence in a stack of N$`\times `$1300 sec exposures we expect 0.02257N pixels at a given position to have suffered a hit. Taking the dispersion to be $`\sqrt{0.02257N}`$ and using $`5\sigma `$ rejection we find that in combining the 98 F606W images the highest 7 pixels should be rejected while for the 148 F814W images we should eliminate the highest 9. Statistically this approach is superior to a median of the frames, which is equivalent to a mean of 2/3 of the number of frames, whereas in our procedure we used 93% of the images in constructing the mean.

For each chip there are 4 combined frames, $`\mathrm{V}_{\mathrm{new}}`$, $`\mathrm{V}_{\mathrm{old}}`$, $`\mathrm{I}_{\mathrm{new}}`$ and $`\mathrm{I}_{\mathrm{old}}`$. A starlist was generated for each chip using DAOFIND on the $`\mathrm{I}_{\mathrm{new}}`$ and $`\mathrm{V}_{\mathrm{new}}`$ frames. We then (1) visually inspected each object in the list and rejected those objects that corresponded to detections of non-astronomical sources (such as diffraction spikes), (2) subtracted a spatially smoothed copy of the image from the image itself and inspected this image for missed stars, and finally (3) combined the 2 lists into our final starlist.

Stars on each of the chips within the 4 frames were photometered using ALLSTAR with the starlist as input. ALLSTAR fits a PSF to each object in the starlist and does recentering so as to achieve the best fit. Note that the ability of ALLSTAR to recenter is essential as this allows shifts between stars on the frames to be measured which is critical when applying the starlist from the new frames onto the first epoch images as the proper motion of the cluster relative to the background amounts to about 1 HST pixel over the 6 year time baseline. To ensure that the individual frames were registered on the cluster stars, through an iterative process, we built a transformation file using only cluster stars (about 240 per chip) which move insignificantly relative to each other as the observed cluster velocity dispersion is 3.5 km/s (Peterson et al. 1996). This amounts to a total internal mean motion of about 0.02 pixels over the 6 years. To compensate for distortions in the HST optics, we used a 20-term transformation equation. This is in lieu of but equivalent to applying the Trauger et al. 1995 corrections. To improve astrometric accuracy, we expanded the pc pixels by a factor of 2 and the wf pixels by a factor of 3 using MONTAGE2. This process is similar to “drizzling” (Fruchter & Hook 2002). The individual frames were then averaged together with high pixel rejection as previously described. It was on these combined, expanded frames that the final photometry was carried out using ALLSTAR. No charge transfer corrections were made as the high background and the uniform distribution of the faint stars across the chips suggested that this was not a serious problem on these images.

The resulting proper motion displacement (PMD) diagram (Figure 1) separates cluster stars from the field population rather cleanly. The clump of non-cluster stars in the PMD diagram clearly has a larger dispersion than that of the cluster; with an internal velocity dispersion of only 0.02 pixels in 6 years, motions within M4 are unresolved to our level of precision. For stars in the V magnitude range 21 to 26, the measured PMD dispersion for the cluster is 0.082 pixels ($`0\stackrel{}{\mathrm{.}}0082`$) which we take as our positional measurement error (confirmed by simulation) while the field clump has a dispersion of 0.277 pixls ($`0\stackrel{}{\mathrm{.}}0277`$). A complete discussion of this population will be given in Rich et al. (2002, in preparation).

In the lower two sections of Figure 1 we display the CMDs for data in different areas of the PMD diagram. The data were transformed to Johnson V and Kron-Cousins I using stars measured in GO 5461 and calibrated with ground-based data as discussed in Richer et al. 1997. This resulted in a significant color term in the transformation of F606W to Johnson V. The transformation equations varied slightly from chip-to-chip but for wf2 it was $`\mathrm{V}=\mathrm{F606W}+0.268(\pm 0.014)(\mathrm{F606W}\mathrm{F814W})+9.401(\pm 0.031)`$. Photometric measurement errors alone for MS stars as returned by ALLSTAR ranged from 0.02 magnitudes at $`\mathrm{V}=25`$ (0.01 at $`\mathrm{I}=22`$) to 0.07 at 29 (0.04 at $`\mathrm{I}=24.5`$). This resulted in $`1\sigma `$ errors in magnitude and color for a MS star at $`\mathrm{V}=29`$, $`\mathrm{V}\mathrm{I}=4.5`$ of 0.10 and 0.12 respectively.

The left panels plot all the objects (stars as well as possible galaxies) measured from all 4 chips with no rejection criteria. The central panels isolate the cluster stars – objects within 0.5 pixel in radius ($`0\stackrel{}{\mathrm{.}}05`$) of the mean motion of the cluster. Our simulations suggest that $`<1`$% of cluster stars will have measured PMDs outside this radius. The right panels contain the inner halo stars, again isolated with a 0.5 pixel radius circle centered on their mean motion with respect to the cluster. Any faint unresolved galaxies will appear in this CMD.

## 3 The Approach to the Hydrogen-Burning Limit

M4 is ideally suited for studies of the pop II hydrogen burning limit, by virtue of its proximity and relatively (amongst globular clusters) low concentration. From the ground, the best attempt to study the lower main sequence was made by Kanatas et al. 1995, who were limited primarily by field star contamination. Deep images using NICMOS and WFPC on HST allowed Pulone, de Marchi & Paresce 1999 to probe estimated masses down to $`0.15\mathrm{M}_{}`$, but with poor statistics. A second epoch of HST observations allowed Bedin et al. 2001 to separate cluster members from background stars using the cluster proper motion, with the faintest cluster members measured at $`\mathrm{F814W}24.5`$ with masses near $`0.1\mathrm{M}_{}`$. Bedin et al. determined a completeness-corrected luminosity function, finding that their last bin at $`\mathrm{I}24.5`$ was statistically indistinguishable from zero. They note that if a larger number of stars are detected and this bin stays empty, then it would strongly indicate the location of the hydrogen burning limit. The considerably deeper observation described above allows us to examine this question in more detail.

Inspection of Figure 1 shows that, while there is indeed a paucity of stars at $`\mathrm{V}>27.5`$, $`\mathrm{I}>23.8`$, there does remain a trail of potential MS stars extending redwards, almost to $`\mathrm{V}=30`$. These stars are present even with very restrictive proper motion cuts and have a high probability of cluster membership (Figure 2) so we are convinced that the LF does not have an abrupt termination bluer than $`\mathrm{V}\mathrm{I}=5`$. Indeed, there is little reason to expect that it will as the current generation of low mass MS models (Baraffe et al. 1997, Cassisi et al. 2000, Montalban et al. 2000) all predict large changes in luminosity and color with very small changes in mass near the hydrogen burning limit. Hence the approach to the hydrogen-burning mass limit should be characterised by an extended LF, such as is indeed seen in M4.

The least luminous MS stars observed in M4 are significantly fainter than any other known subdwarf at this metallicity (Leggett et al. 1998a; compare with LHS 1742a and LHS 377). An immediate consequence is that there should be correspondingly faint red subdwarfs in the field waiting to be discovered, although they will be rare. For every field subdwarf in the range $`7.5<\mathrm{M}_\mathrm{V}<15`$, there should be $`0.02`$ subdwarfs with $`15<\mathrm{M}_\mathrm{V}<17.4`$, or less than 1 in 30,000 of all stars in the Solar neighbourhood.

There is also a small sample of apparent cluster members located between the cooling sequence and the MS, and some are fainter in V than even the least luminous WD detected. It is possible that these are field objects with proper motions similar to the cluster. If they are cluster members, they could be low luminosity He WDs in binaries with extremely low mass MS stars. Indeed, such objects are predicted to be the remnants of cataclysmic variables currently accreting at very low rates (Townsley & Bildsten 2002). Alternatively, these may indicate a range of colors for pop II stars near the hydrogen-burning limit, since low mass, low metallicity stars can be quite blue (Saumon et al. 1994).

## 4 The Cluster Mass Function

The derivation of a MF from the LF for M4 is an uncertain procedure because the current theoretical models for low temperature stars of this metallicity do not fit the lower MS very well. Modelling difficulties are due to problems in treating convection and opacity in the optically thin, dense and cool outer atmospheres of these stars where molecules dominate and where non-ideal gas effects are important. This affects both the color-Teff relations and the stellar luminosity. At the lowest masses severe problems with calculating the adiabatic gradient occur with resultant uncertainty in the mass-luminosity relation (Montalban et al. 1999). With these caveats, we generated and display the cluster MF in Figure 3 using the models of Montalban et al. to convert luminosities to masses. The input luminosity function was corrected for incompleteness by inserting in 4 trials a total of 4900 stars with magnitudes from $`\mathrm{V}=20`$ to 31 to each chip. For an artificial star to be considered recovered it had to be found on at least 3 of the frames (so that its motion could potentially be determined) as well as being recovered within 1 magnitude of its input brightness and within 0.5 pixels of its input position.

The resulting MF rises very slowly to low masses, with a formal slope $`\alpha =0.75`$ ($`\mathrm{N}(\mathrm{m})\mathrm{m}^\alpha `$, Salpeter $`\alpha =2.35`$). A $`\chi ^2`$ fit finds that slopes of $`\alpha =0.55`$ – 1.05 are acceptable at the 99% confidence level, a tighter constraint will be possible when better lower MS models are available. There is no evidence that the MF turns over, even to the lowest masses we detect ($`0.09\mathrm{M}_{}`$). This assertion is more solid than might appear from the uncertain model fitting. The total number of stars fainter than $`\mathrm{I}21`$ (brighter than this the models do fit the data) is consistent with the expectations based on the extrapolated MF down to $`\mathrm{M}0.09\mathrm{M}_{}`$, regardless of the M-L relation. We can also predict the number of MS stars expected to $`0.08\mathrm{M}_{}`$ (the brown dwarf limit for solar metallicity). Extrapolating the current MF suggests that we should find 20 more stars redder than $`\mathrm{V}\mathrm{I}=5`$. If the hydrogen-burning limit for M4 is at higher mass, e.g. at $`0.086\mathrm{M}_{}`$ as suggested by the Montalban et al. models, this will drop to 8. Within this context it is of interest to note that we found no stars with I magnitudes in the range 25 – 27 that possess PMDs appropriate for cluster membership and are found solely on the I frames.

The WD population also allows us to constrain the MF slope above the cluster turnoff as discussed in Richer et al. 1997. The completeness-corrected WD counts are 602 to $`\mathrm{V}=30`$. The detected 570 MS stars between $`0.09\mathrm{M}_{}`$ and $`0.65\mathrm{M}_{}`$ (upper limit due to saturation) imply 1219 WDs if the slope $`\alpha =0.75`$ extends all the way through the turnoff mass ($`0.8\mathrm{M}_{}`$) to the expected upper limit for WD progenitor masses ($`8\mathrm{M}_{}`$). Thus, we find a factor of two less than expected. However, if the MF slope is at the steep end of the acceptable range found above, then we predict 584 WDs, in agreement with the observed numbers. However, the WD counts are far above that expected if the global M4 IMF of Pulone et al. 1999 is adopted. They propose $`\alpha 2.4`$ above the turnoff mass which would predict a paltry 264 WDs. Furthermore, if M4 has undergone significant mass segregation as they suggest, this number would be even smaller.

If there is a significant fraction of the cluster WD population which possess pure He atmospheres, they will have cooled beyond our detection limit (Hansen 1999) and will not be included in the above census. A population comprising 50% of all cluster WDs would allow a continuous power law $`\alpha =0.75`$ from $`0.09\mathrm{M}_{}`$ to $`8\mathrm{M}_{}`$. Hence the constraints from WD counts that can be placed on the MS are dependent upon assumptions about the shape of the MF and the WD fraction with He atmospheres. However, the fact that all of the observed WDs seen can be accounted for with a MF slope of $`\alpha 1`$ is at least consistent with a simple picture where most WDs are H-rich, and where we have indeed seen the termination of the WD cooling sequence in M4.

HBR and GGF are supported in part by the Natural Sciences and Engineering Research Council of Canada. HBR extends his appreciation to the Killam Foundation and the Canada Council for the award of a Canada Council Killam Fellowship. RMR, MS and IS acknowledge support from proposal GO-8679 and BH from a Hubble Fellowship HF-01120.01 both of which were provided by NASA through a grant from the Space Telescope Science Institute which is operated by AURA under NASA contract NAS5-26555. BKG acknowledges the support of the Australian Research Council through its Large Research Grant Program A00105171.