text:
# 1 Flow Chart for Sequential Data Assimilation

A Kalman Filter for Ocean Monitoring

Konstantin P.Belyaev

LNCC, Petropolis, Brazil

Detlev Müller

MPIMET, Hamburg, Germany

Abstract: The feasibility of global ocean state estimation by sequential data assimilation is demonstrated. The model component of the assimilator is the GROB version of the MPIMET ocean circulation model HOPE. Assimilation uses the Fokker-Planck representation of the Kalman Filter. This approach determines the temporal evolution of error statistics by integration of the Fokker-Planck Equation. Phase space advection and diffusion are obtained from histogram techniques considering the model as a black box. For efficiency, the estimation procedure utilizes a combination of nudging and Kalman Filtering. The ocean state is estimated for the El Nino year 1997 by dynamical extrapolation of observed sea-surface temperatures and TAO/TRITON subsurface temperatures. The model-data combination yields improved estimates of the ocean’s mean state and a realistic record of El Nino related variability. The assimilator proves as an efficient, viable and thus practical approach to operational global ocean state estimation.

1.) Introduction

Ocean state estimation draws its wider societal as well as scientific significance from the central role of the oceans in Earth’s climate system. At a time when the impact of climate variability on societal infrastructures is increasingly felt, the need for comprehensive climate monitoring is generally accepted. While mankind is primarily affected by meteorological manifestations of climate variability, large-amplitude weather fluctuations oftentimes screen the atmospheric climate signal. In practice, atmospheric climate observation proves prohibitively intricate. Alternatively, estimates of the state of the ocean interior with its enormous capacity to store and distribute water, heat and radiatively active trace substances (such as carbon dioxide) provide direct evidence of the climate signal. As the dominating climate component, the ocean acts as a Markov Integrator of atmospheric noise, provides the memory of the climate system and sets time-scales of climate processes by (at least partly) predictable transport mechanisms. Thus, practical climate monitoring anchors on operational ocean state estimation.

Only a decade ago, an observational basis for the assessment of a temporally changing global ocean circulation was practically nonexistent. In fact, oceanography was plagued by a notorious undersampling problem and the ocean circulation was widely considered as a stationary flow. The most comprehensive data set of this era is the World Ocean Atlas (WOA) by Levitus and coworkers , a meticulous collection of hydrographic data from oceanographic archives all over the world. This atlas represents a global ocean density field that is best associated with the mean ocean circulation during the second half of the twentieth century. Temporal changes such as those at the roots of past, present and future climate variability are far beyond the scope of these data. Since then, ocean observation underwent an explosive development. Today, space-borne ocean observatories provide global data sets of the mesoscale state of the sea surface in near real-time. Oceanographic data paucity has been replaced by an almost overwhelming data stream . Nevertheless, observational monitoring of the ocean interior remains technically difficult and costly for the foreseeable future. For mere book keeping as well as for analysis and interpretation, ocean state estimation thus relies heavily on the dynamical extrapolation of observations by means of numerical models of the global ocean circulation.

Such models play a central role in the study of the climate system and its ocean component for more than three decades. As engineering devices with the capacity for otherwise impossible experimentation they provide a laboratory for the simulation, extrapolation and ultimately understanding of the observational record. So far, the most convincing demonstration of the potential of circulation modeling has been the discovery of the El Nino mechanism and its repercussions on remote regions of the globe. After Matsuno’s theoretical ground work  and the provision of observational evidence by Wyrtki , the essential dynamical processes were numerically simulated with simple shallow water models . Today, operational El Nino monitoring and forecasting are implemented at a number of institutions world wide and the variability of tropical circulations has become one of the best understood aspects of the climate system.

On the other hand, dedication of more complex models to the same problem has yet been unable to advance El Nino simulation and forecasting significantly. Nor has it been possible to gain comparable insights for the extratropics. While scientists are well aware of extratropical circulation variability such as the North Atlantic or Arctic Oscillation or the Antarctic Circumpolar Wave, numerical analyses have clearly been less yielding in these instances. Uncertainties of this kind also raise the question of the adequacy of contemporary numerical circulation models for longer term climate projections.

The major difficulty of global circulation modeling is the lack of a physically consistent and numerically soluble formulation of the circulation problem. Physically, the global ocean circulation poses the thermohydrodynamical problem for a viscous fluid on the rotating geoid. The gravest problem of circulation theory is certainly the absence of an energetically consistent formulation of the equations of motion for a viscous rotating fluid. While such formulations are well known for viscous nonrotating fluids and ideal rotating fluids, the derivation of the energy budget from the equations of motion of a dissipative rotating fluid remains a matter of scientific debate. As one consequence, wave-dissipation in rotating fluids is not very well understood. Moreover, viscosity plays a significant role for the numerical stability of circulation models. It is clearly desirable that parameterizations of subscale transports steer momentum, energy and vorticity along realistic, i.e. energetically and vortically consistent pathways in space-time and wave vector space. For contemporary circulation models, this matter remains essentially unsolved.

Satellite orbits reveal a complex fine structure of the planet’s shape and considerable uncertainties particularly about the marine geoid still exist. For most purposes of circulation modeling, however, an approximation in terms of a spheroid or even a sphere is probably sufficient. A source of dissatisfaction with 3-dimensional spheroidal or spherical circulation equations are difficulties in analytically obtaining simple stationary solutions which reflect characteristic circulation features such as geostrophy and thermal wind balance. Moreover, linearizations in these coordinates do generally not admit separation of variables. Hence, it has yet been impossible to study analytically the propagation of acoustic, gravity and Rossby wave disturbances together with the stability of simple flows in a common, 3-dimensional framework. For numerical integration, such problems are insignificant. Nevertheless, most contemporary numerical circulation models compromise geometric-dynamic integrity in favor of a multi-$`\beta `$-plane approximation to the geometry of the geoid. This approach codes Laplace operators as sum of second order derivatives and ignores first order contributions from nontrivial Christoffel symbols.

A third set of issues of circulation modeling is associated with the nonlinearity of fluid dynamics. Nonlinear field theories inevitably couple variability on the smallest space- and time-scales with the largest scales available. For numerical as well as theoretical purposes, processes on small spatial and fast temporal scales should be eliminated from circulation equations. In the first place, this applies to acoustics. The widely favored sound filter invokes the Boussinesq approximation to inertia and weight of the fluid. While this approach has been quite successful in the study of internal gravity waves and Rayleigh-Benard instability, its vorticity-inconsistency becomes a problem in long-term integrations of circulation modeling. Furthermore, Rayleigh-Benard (or static) instability involves fast convective motions on small spatial scales. Such convective events are generally not resolved in global circulation models and their crucial role for the thermohaline circulation is represented by appropriate parameterizations. Hence, models assume the ocean to be in hydrostatic equilibrium and simply neglect internal vertical accelerations relative to Earth’s gravitational acceleration. However, this straightforward introduction of hydrostatics generates a problem. Now, momentum density is a 2-dimensional vector while the mass flux vector of the continuity equation remains 3-dimensional. Such a violation of the first law of motion will generally be uncritical for stationary flows or order-of-magnitude estimates. However, for the long-term integrations of circulation modeling, such a formulation is not entirely satisfactory.

The Boussinesq approximation, hydrostatics and the so-called “traditional approximation” of inertial forces define Richardson’s Primitive Equations . Currently, these equations provide the physical basis for most global circulation models. While the Primitive Equations account consistently for velocities and thermodynamics of equilibrium circulations, they do not pose a Newtonian dynamical problem.

The present study utilizes GROB HOPE, a coarse version of the numerical Hamburg Ocean Primitive Equations model of MPIMET . This model is based on a C-grid discretization of the Primitive Equations for UNESCO seawater  and allows various convection parameterizations which account for different characteristics of this process in the open ocean and in the bottom boundary layer down submarine slopes. GROB HOPE includes sea ice dynamics with viscous-plastic rheology  parameterizing cracking, ridging, rafting and deformation of sea ice. The model is forced by buoyancy fluxes and wind stresses at the sea surface as well as the freshwater discharges of Earth’s 50 largest rivers.

In long term experiments (integration time: 1000 years) with climatological forcing resolving the annual cycle, the model assumes an essentially drift-free cyclostationary state after a few centuries which reproduces the major water masses and gyre structures of the global ocean circulation as well as the sea ice cover and its seasonal variation at high latitudes. While this model circulation exhibits the characteristic degree of realism of state-of-the-art simulations it also displays a number of typical deficits. The model fails to maintain the observed Pacific Intermediate Waters. Furthermore, while the poleward Atlantic heat transport is certainly of the observed order of magnitude, its maximum of 0.8 PW is still somewhat lower than the 1.1 PW suggested by observations. On the other hand, the mass transport by the Antarctic Circumpolar Current with 180 Sverdrup in the Drake Passage is higher than the observed 140 Sverdrup. The path of the Gulf Stream which is crucial for the European climate and weather turns out to be quite sensitive to the details of the atmospheric forcing and the chosen parameterization of subscale transports. For an extensive discussion of the strengths and weaknesses of the GROB HOPE circulation, see . At this time, the versatility of numerical circulation models is sufficiently developed to simulate a wide variety of preconceived scenarios. However, beyond this illustrative role, state-of-the-art models generally lack the capacity of scientific discrimination between competing predictions and projections. To a large part, these uncertainties can be resolved by the systematic combination of models with observations.

Such combination and confrontation of models with extensive and novel observational data sets has been the outstanding factor in forecast-improvement by numerical weather prediction over the last thirty years. With maturing numerical ocean circulation models and a growing understanding of the variability of the ocean circulation this approach has now also become attractive to the oceanographic and climate communities. Given the characteristics of model data and observations in Earth System Modeling, the integration of information from different sources poses a considerable data engineering problem. Frequently, observations do not refer to prognostic model variables while other parameters such as vertical velocities are practically unobservable. Moreover, observations are typically distributed highly irregular in space-time. And data sets from both sources, model and observation, are large. The mathematics of the optimizing synthesis of large data sets are the objective of estimation theory . To avoid the mutual enhancement of model- and data-error, estimation theory has been (and still is) developing a number of what are called data assimilation algorithms. Generally, these algorithms fall into two classes: variational and sequential techniques . The equivalence of both methods is readily demonstrated in simple cases.

Variational assimilation, namely the Adjoint Method, is based on an application of inverse modeling techniques to the estimation problem. Variation of control parameters minimizes a cost function formed by the model-data misfit. This approach lends itself particularly to the estimation of equilibrium states and processes of finite duration. Computation of the cost gradient with respect to the controls calls for what is often referred to as the temporally backward integration of the adjoint model. For complex models, coding of the model adjoint is a substantial task, well comparable to coding the model itself. The practical relevance of adjoint assimilation in Earth System Modeling therefore arose only after the advent of the theory of automatic differentiation  and the subsequent development of automatic adjoint code compilers . With global circulation models, the Adjoint Method finds presently wide application in sensitivity studies.

Sequential methods such as the Kalman Filter are more specifically tailored to the needs of monitoring and prediction . These updating schemes emerge from the application of the theory of stochastic processes to estimation and yield an estimate of minimum variance. On update, the relative weight of model and data is determined by the Kalman Gain which is computed from data and model error dynamics. To this end, the model error is considered as a stochastic process. The dynamics of such processes can be equivalently formulated in the Langevin (or Heisenberg) representation and in the Fokker-Planck (or Schrödinger) representation . The Langevin picture addresses the space-time behavior of the process in terms of its moments. In practice, this refers generally to the covariance only. Formally, the temporal development of the covariance is uniquely determined by the model dynamics. However, the practical derivation of the covariance dynamics for a complex model such as a global ocean circulation model readily becomes everything but straightforward. This applies particularly to nonlinear models, the issue of boundary and initial conditions for the covariance, stability questions and the problem of temporally backward assimilation. Nevertheless, at this time the literature on Kalman Filter assimilation in Earth System Modeling and other branches of engineering is almost exclusively dominated by the Langevin approach .

Alternatively, a stochastic process may be considered in phase-space in terms of its probability density. Provided the process is Markovian and jumps remain small in an appropriate sense  the dynamics of this probability density are governed by the Fokker-Planck Equation. The advection- and diffusion-coefficients of this linear parabolic differential equation are determined by model dynamics and observational error statistics. In general, these coefficients are also difficult to obtain from a complex model. However, for sufficiently short update intervals, phase-space advection and diffusion can be determined phenomenologically from the model output by histogram techniques. In this framework, the assimilation method provides practical answers to the issues of phase-space reduction, model nonlinearity, initial and boundary conditions for higher moment dynamics as well as stability. Moreover, the existence of the Backward Fokker-Planck Equation  will permit the generalization of sequential assimilation to include the temporally backward extrapolation of data information. The mathematical aspects of the Fokker-Planck representation of sequential Kalman Filter assimilation have been developed in detail by .

With the typical volume of model output and observational record in Earth System Modeling, computational demands for assimilation with least-square optimality are always quite high. For a reduction of the computational burden, the present estimation utilizes a combination of Kalman Filter assimilation and simple “nudging”. While subsurface temperatures from the TAO/TRITON array will be assimilated sequentially, observations of global sea-surface temperatures are essentially inserted into the model at daily intervals. The feasibility of this simplistic technique is by no means trivial. Older model generations were generally unable to “digest” essentially unprocessed data and model-data inconsistencies would readily emerge in various regions of space-time and phase space. It will here be shown that the quality of contemporary models and data sets is sufficiently high for nudging to be beneficial for the ocean state estimate.

2.) Fokker-Planck Picture of the Kalman Filter.

Societal and scientific needs in Earth System Monitoring are presently best met by an operational, steady combination and confrontation of substantial, yet incomplete observations with complex, but nevertheless approximate numerical models. This combination aims to fill data gaps by dynamical extrapolation of observations on the basis of physical laws incorporated in the model and simultaneously constrain model uncertainties by operating the model in close vicinity of observations. One source of information is used to compensate the deficits of the other and thus arrive at a comprehensive state estimate including an assessment of model- and observation-quality. The assimilation algorithms of estimation theory are information integration techniques which prevent the mutual enhancement of errors from different sources and maximize the benefit from imperfect models and data. As an engineering tool, data assimilation accepts or rejects a hypothesis (the model) while a constructive evaluation of structural model deficits remains beyond its scope.

The key idea of sequential assimilation is to integrate the model until an observation becomes available. At this time, model integration is halted and the state of the system is updated by an appropriate combination of model prediction and observation. Subsequently, this update provides the initial value for the continued model integration (fig.1). Hence, the main task in sequential data assimilation is the determination of the temporal development of the relative weight of model prediction and observation.

The state of the system under consideration (here: the global ocean) is given by a finite- (generally: very high) dimensional vector $`Y(t,x)`$ where $`t`$ denotes time and $`x=(x_{(1)},x_{(2)},x_{(3)})`$ the 3-dimensional spatial coordinate. In numerical ocean models, time and space coordinates will generally be members of a discrete, 4-dimensional grid. The time-development of the state of the system is governed by the model dynamics

$$\frac{d}{dt}Y_m(t)=\mathrm{\Lambda }(Y_m,t)$$
(1)
where the index m refers to the model and $`\mathrm{\Lambda }`$ is a generally nonlinear operator of $`Y_m`$ with

$$\mathrm{\Lambda }(Y_m=0,t)=0$$

and some explicit time-dependence representing external forcing. For simplicity, spatial coordinates have been suppressed in (1). It is now assumed that the model is sufficiently “good” so that the time-development of the true, i.e. observed system is given by the equations of motion

$$\frac{d}{dt}Y_o(t)=\mathrm{\Lambda }(Y_o,t)+W_0$$
(2)
where the index $`o`$ refers to observation, while $`\mathrm{\Lambda }`$ is the same operator as in (1) and $`W_0`$ an additive noise of known probability distribution. The possibly nonstationary noise is assumed to have zero average, finite spatial and no temporal correlation, i.e. it is assumed to be temporally white

$$<W_0(t,x)>=0,<W_0(t,x)W_0(t+\tau ,x+r)>=Q(t,r)\delta (\tau ).$$

Physically, this noise accounts for model deficits. Equation (2) is the Langevin representation of a (nonlinear) stochastic differential equation.

The model error

$$y(t,x)=Y_o(t,x)Y_m(t,x)$$

appears under a variety of names in different applications of assimilation theory. Variational techniques generally term this quantity “misfit” while numerical weather prediction refers to the model error as “innovation”. By definition, the model error satisfies the Langevin equation

$$\frac{d}{dt}y(t)=\mathrm{\Lambda }(y+Y_m,t)\mathrm{\Lambda }(Y_m,t)+W_0.$$

With the help of the mean value theorem, the dynamical operator is rewritten as

$$\mathrm{\Lambda }(y+Y_m,t)\mathrm{\Lambda }(Y_m,t)=y\mathrm{\Lambda }^{}(Y_m,t)+R$$

where the residual R comprises terms higher than second order. Writing now

$$\mathrm{\Lambda }(y+Y_m,t)\mathrm{\Lambda }(Y_m,t)=\mathrm{\Lambda }(y,t)\mathrm{\Lambda }(y,t)+y\mathrm{\Lambda }^{}(Y_m,t)+R$$

and using again the mean value theorem on $`\mathrm{\Lambda }(y,t)`$ together with $`\mathrm{\Lambda }(0,t)=0`$, one arrives at

$$\frac{d}{dt}y(t)=\mathrm{\Lambda }(y,t)+W$$
(3)
where the new noise $`W`$ comprises $`W_0`$ and terms of second and higher order. It is readily seen that

$$<W(t,x)>=0.$$

The error dynamics (3) account for nonlinearities and provide the basis for the histogram technique to be invoked below.

For simplicity, data are here assumed to represent prognostic model variables directly. Furthermore, measurements are made at specific observation points $`X=(X_{(1)},X_{(2)},X_{(3)})`$ which are generally much less in number than model grid points $`x`$ and do not coincide with these. Interpolation to model points will generally be straightforward and observation points are here assumed to coincide with model points. In the context of monitoring and prediction, observations typically become available at specific times which are termed “analysis times” in numerical weather prediction. At these times, sequential data assimilation halts the model integration and updates the system’s state vector with observational information. This update is given in terms of a convolution-type integral

$$Y(t,x)=Y_m(t,x)+_0^t𝑑\tau 𝑑XG(\tau ,x,X)y(\tau ,X)$$
(4)
where the kernel $`G(\tau ,x,X)`$ is called the gain. For the update (4), the gain accounts for all observations prior to analysis time and hence only represents the temporally forward propagation of data information. In this sense, the kernel in (4) represents the “retarded gain”. It is equally possible to ask: at which state did the system initially start, if it is at the observed state on analysis time. The account of the temporally backward propagation of data information requires the inclusion of the “advanced gain”. Sequential assimilation admits the dynamical extrapolation of data information into both, past and future. The corresponding algorithm, namely a meaningful representation of the advanced gain becomes particularly transparent in the phase space representation. Its formal details will be discussed elsewhere.

There is a number of expressions for the gain-kernel in (4) which all lead to a different variance for the state estimate. The estimate has minimum variance if the gain is given by the Wiener-Hopf equation

$$K(t,x,X)=_0^t𝑑\tau 𝑑X_0G(\tau ,x,X_0)K(t\tau ,X,X_0)$$
(5)
where

$$K(t,x_1,x_2)=<y(t,x_1)y(t,x_2)><y(t,x_1)><y(t,x_2)>.$$

is the error covariance with brackets indicating the ensemble mean and unbraced indices denoting different points, i.e. different 3-tuples (and not different vector components). The gain defined by equation (5) is called the Kalman Gain and the sequential data assimilation algorithm given by (4) and (5) is called the Kalman Filter . It follows from (5) that the central problem of Kalman Filter assimilation is the determination of the temporal development of the error covariance $`K(t,x_1,x_2)`$.

Formally, the covariance matrix and its temporal behavior are uniquely determined by the model error dynamics (3). In practice, however, exploitation of this equation encounters a number of serious problems. Covariance dynamics on the basis of (3) are determined at every model grid point. If model integration requires $`N`$ operations per time-step, integration of the covariance dynamics requires additional $`N^2`$ operations. An order of magnitude estimate of the number of operations for a contemporary ocean circulation model such as GROB HOPE is $`N10^5`$. Even taking into account the rate of increase in available computing capacities, integration of the complete covariance dynamics in the Langevin picture becomes prohibitive for global circulation models. Furthermore, these models typically involve strong nonlinearities. For mean quantities such as the the error covariance, this nonlinearity leads to a hierarchy problem since the mean of a product does not equal the product of the means. The difficulties and ambiguities of practical closures of such hierarchies are well known from turbulence theory . Additional complications arise from the lack of initial and boundary conditions for the dynamics of higher moments and for the determination of stability and uniqueness conditions of the assimilation procedure. These problems represent some of the reasons for the prevalence of the Adjoint Method in Earth System estimation.

An alternative approach to the Kalman Filter utilizes the phase space or Fokker-Planck representation of stochastic processes. In this framework, a detailed assimilation algorithm for application in Earth System estimation has been developed by Belyaev and coworkers . The starting point of this formulation is the joint probability distribution

$$p=p(t,\eta _1,\eta _2;X_1,X_2)$$

for the two-component projection $`\eta `$ of the model error $`y`$

$$\eta =(\eta _1,\eta _2)=(y(t,X_1),y(t,X_2))$$

to have the value $`\eta _1`$ at observation point $`X_1`$ and the value $`\eta _2`$ at observation point $`X_2`$ at time t. In terms of this probability the error covariance (5) takes the form

$$K(t,X_1,X_2)=d^2\eta \eta _1\eta _2pd^2\eta \eta _1pd^2\eta \eta _2p.$$
(6)
Determination of the error covariance thus requires the calculation of joint probability distributions p and this operation sets the magnitude of necessary computing capacities. Notice that (6) defines the error covariance only for pairs of observation points. In Earth System estimation, the number $`M`$ of these points is typically orders of magnitude smaller than the number $`N`$ of model grid points. Using the symmetry of the joint probability distribution, the covariance matrix for $`M`$ observations has $`M(M+1)/2`$ independent members. Given 1000 data points (which still is a fairly small data set in Earth System observation), the error covariance calls for the computation of half a million probability functions. While this number is clearly much smaller than $`N^210^{10}`$, it still poses a considerable computational task.

In the phase space picture, the error $`\eta `$ is viewed as a Markov process. The temporal development of the conditional probability distribution of such a process is governed by the Master Equation . For Markov processes with small jumps this equation is well approximated by the Fokker-Planck Equation . In the present case this equation takes the form

$$_tp=^n\mathrm{\Lambda }_np+\frac{1}{2}^n(q_{mn}^mp)$$
(7)
where the usual summation convention for indices $`m,n,\mathrm{}=1,2`$ is implied. The vector

$$J_n=\mathrm{\Lambda }_np\frac{1}{2}q_{mn}^mp$$

denotes the advective-diffusive probability flux in phase space with advection velocity $`\mathrm{\Lambda }_n`$ and diffusion $`j_n=\frac{1}{2}q_{mn}^mp`$. In this sense, the Fokker-Planck Equation (7) expresses the conservation of probability in phase space. Temporally, this linear parabolic equation governs the development of the conditional probability $`p(t_2,\eta _2|t_1,\eta _1)`$ of the error $`\eta `$ to have the value $`y_2`$ at time $`t_2`$ and location $`X_2`$ given its value $`y_1`$ at time $`t_1<t_2`$ at location $`X_1`$. The required joint probability follows from the solution of the Fokker-Planck Equation by Bayes’ Rule

$$p(t,\eta _2,\eta _1)=p(t,\eta _2|t,\eta _1)p(t,\eta _1).$$

It is noted that there is also a Backward Master Equation . The small jump approximation to this equation leads to the Adjoint Fokker-Planck Equation which governs the temporally reversed development of the error distribution. Hence, the Adjoint Fokker-Planck Equation provides the basis for the determination of the Advanced Kalman Gain.

For the solution of the Fokker-Planck Equation (7) the phase space advection $`\mathrm{\Lambda }_n(t,\eta )`$ and the diffusion tensor $`q_{mn}(t,\eta )`$ have to be known. These parameters are determined by error dynamics and data stochastics according to

$$\mathrm{\Lambda }_n(t,\eta )=<\mathrm{\Lambda }_n(t,\eta )|\eta =y>=\tau ^1d^2\eta ^{}(\eta _n\eta _n^{})p(t,\eta |t^{},\eta ^{})$$

with $`\tau =tt^{}`$ for $`t^{}<t`$ for the advection while one has for the diffusion tensor

$$q_{mn}(t,\eta )Q_{mn}(t,\eta )=<\mathrm{\Lambda }_m(t,\eta )\mathrm{\Lambda }_n(t,\eta )|\eta =y>=\tau ^1d^2\eta ^{}(\eta _m\eta _m^{})(\eta _n\eta _n^{})p(t,\eta |t^{},\eta ^{})$$

where $`Q_{mn}`$ denotes the data covariance. With these definitions, the Fokker-Planck Equation (7) is readily seen to be equivalent to the Langevin Equation (3). Multiplying (7) from the left by $`\eta `$ and integrating over phase space one obtains

$$\frac{d}{dt}<\eta _n>=<\mathrm{\Lambda }_n(t,\eta )>$$

in agreement with the two-component projection of the ensemble average of (3).

In principle, the advection- and diffusion-terms are determined by the Langevin representation (3) of the error dynamics. However, particularly for strongly nonlinear dynamics such a derivation encounters serious formal difficulties and a unique and practical solution to this problem is currently not known. In view of the typical problem in Earth System Monitoring, Belyaev and coworkers  propose an alternative, phenomenological determination of these parameters. In Earth System Monitoring the time-interval between consecutive samples is typically short compared to the time scales of the processes under observation. Under these conditions it becomes possible to consider the model as a black box and determine the transition (i.e. conditional) probabilities by histogram techniques from model input and output. To this end, the number $`N^{}`$ of all grid points is counted for which $`\eta `$ has the value $`y^{}`$ at time $`t^{}`$. At the later time $`t=t^{}+\tau `$ all former $`y^{}`$-points $`N`$ are counted for which $`\eta `$ now has the value $`y`$. The conditional probability is then given by the ratio

$$p(t,\eta |t^{},\eta ^{})=N/N^{}.$$

Since $`0NN^{}`$, this expression always satisfies the condition

$$0p(t,\eta |t^{},\eta ^{})1$$

necessary for being interpreted as a probability density. With the help of this probability, the advection and diffusion parameters are readily obtained by phase space integration according to the above formulas.

It is now possible to solve the Fokker-Planck Equation numerically. As a linear parabolic differential equation in the 2-dimensional unbounded plane, the equation is efficiently integrated by the Peaceman-Rachford scheme . For details of this integration concerning initial conditions, positive definiteness and normalizability of the solution see . The resulting probability density determines the error covariance at all observational points $`X`$ according to (6) and the covariance at all model grid points is constructed from this expression by interpolation . Using this error covariance, the Wiener-Hopf Equation (5) is solved for the Kalman Gain and the model is finally updated according to (4).

3.) Ocean State Estimation

The feasibility of operational global ocean state estimation will here be demonstrated by combining simulations of the numerical circulation model GROB-HOPE with observations of global sea-surface temperatures (SST) and observed subsurface temperatures from the TAO/TRITON array for the El Nino year 1997. Besides a globally realistic mean state, the objective of the estimate is the improvement of the model’s El Nino simulation.

GROB HOPE has 20 layers in the vertical with high resolution of 10 layers in the upper 500m. In the horizontal, the model uses a spatially inhomogeneous grid obtained from a conformal transformation of the geographical coordinates. At the present stage of model development and availability of computer capacities, the disadvantages of an inhomogeneous horizontal grid are easily outweighed by its advantages. For one, polar singularities are avoided by transformation of the model poles to a continental site. Secondly, the spatial inhomogeneity of the horizontal grid allows high resolution in regions of interest (up to 25 km for the Arctic Ocean in the present case) while low resolution is accepted for remote regions (300 km near the equator in the present case). This design avoids well-known open boundary problems of fine-resolution regional or nested models. While the low-resolution regions provide a model-consistent climatology, the high-resolution regions admit even the study of mesoscale processes. In spite of this versatility, the machine requirements for GROB HOPE are those of a global model with a spatially homogeneous $`3^{}\times 3^{}`$ grid. This design permits a time step of 2.4 hours. With its coarse spatial resolution in the tropics, the GROB version of HOPE does not especially qualify for El Nino simulation. It is here to be shown that assimilation of observations is able to offset these design limitations. Success in this framework provides a demonstration of the capacities of sequential assimilation. For operational purposes, on the other hand, data will always be assimilated into the best model available.

After an initial spin-up period of 2 years with restoring to the 3-dimensional buoyancy climatology of WOA the model is integrated from 1948 to the present with surface forcing derived from the NCEP reanalysis . Atmospheric data are interpolated onto the GROB HOPE grid and surface buoyancy- and momentum-fluxes are calculated by bulk formulae  depending on both, the atmosphere and the ocean. Hence, the eventual ocean forcing is determined by the particular realization of the ocean state by the model while the present ocean-only set-up is unable to account for a feedback of the ocean on the atmosphere. For a reduction of trends in the deep

ocean, integration over the NCEP period is repeated. Furthermore, model surface salinities are nudged to a mean annual cycle taken from WOA with a time constant of a little over a year (385d). Use of a mean annual cycle rather than an annual mean accounts for the seasonal variation in the hemispheric distribution of convective activity. With this forcing the model is integrated to 31 December 1997. The period from 1 January 1997 to 31 December 1997 is taken as the control run in the present experiment and the model state at 31 December 1996 provides the initial condition for the assimilation. It is noted that model runs considered here do not address the prediction problem. Surface data transfer external El Nino information to the ocean model.

Fig.2 shows the monthly mean of the net surface heat flux of the control configuration for December 1997. This heat flux is determined by atmospheric data from the NCEP reanalysis and oceanic data from GROB HOPE. The main feature is the characteristic seasonal separation of the (southern) summer- and (northern) winter-hemisphere: the ocean gains heat in summer and loses heat during winter. A particular detail in the North Atlantic is associated with model problems in simulating a realistic Gulf Stream path: off the American east coast, the ocean is unrealistically warm leading to a pronounced heat loss while the ocean is unrealistically cold in the region of the so-called North West Corner leading in turn to a pronounced heat gain by the ocean. Similar aberrations are seen in the Kuroshio region, the confluence of the Malvinas and Brazil Currents off the South American east coast and for the Agulhas Current near the Cape of Good Hope. The paths of these currents are essentially determined by vorticity dynamics and mismatches of NCEP derived forcing and model simulation are due to ambiguities in the vorticity dynamics of the Primitive Equations. Given the NECP fluxes, GROB HOPE fails to simulate mesoscale details of the state of the underlying ocean surface realistically.

It will now be shown that nudging of observed SST into the model improves the state estimate considerably. To this end, GROB HOPE is restarted from 31 December 1996 and daily mean Reynolds SST of the NCEP data set are inserted into the model’s top layer at a time constant of one day. During this one year integration, model-data incompatibilities do not develop. This is also true for GROB HOPE runs with SST nudging over the full NCEP period (not shown). Fig.3 depicts the monthly mean net surface heat flux for December 1997 with SST nudging. In comparison to fig.2, it is seen that aberrations of the major current systems are significantly reduced and the estimate of mesoscale features of the state of the sea surface improves without penalty.

Nudging effects are not confined to the upper ocean alone. In convectively active regions surface temperature information is rapidly communicated to the abyss. For the present integration period of one year the deep ocean remains of course unable to adjust to the “injected” information. Nevertheless, with these data and for this model, nudging becomes a practical option of ocean state estimation by an efficient and yet robust model-data combination. Other presently available observations are of similar quality: sea-level data from space-borne altimeters and space-based observations of sea-ice cover. By nudging observations of this type into a global ocean circulation model, it is currently possible to arrive efficiently at a comprehensive and realistic estimate of the global state of the sea surface at mesoscale resolution.

For the assessment of the state of the interior ocean consider the equatorial temperature field during the El Nino episode of 1997/98. Fig.4 shows the temperature difference Nudge-WOA along the equator for December 1997 where “Nudge” refers here to the GROB HOPE simulation of the global ocean circulation with NCEP forcing and nudging of daily SST observations, i.e. the run also portrayed in fig.3.

In the abyssal Pacific, simulation and observation are seen to differ by typically less than half a degree. While the simulation is systematically colder than WOA, structural mismatches do not emerge. The agreement is less satisfactory in the abyssal Indic and Atlantic. In the near-surface Pacific, the model clearly exhibits the characteristic El Nino pattern. Relative to the WOA climatology, the eastern and central Pacific are colder while the West is anomalously warm. Comparison with observed subsurface temperatures  shows that the model simulates the phase of the process quite realistically. Since phase information is directly provided by forcing data, this model response is primarily indicative of the consistency of the simulation of near-surface wave propagation with with surface boundary conditions.

Other features of fig.4 exhibit a lesser degree of realism. The warm anomaly in the surface waters of the central Pacific cannot be found in the observational record . Here, the mixed-layer model of GROB HOPE fails to mix the heat supplied at the surface, sufficiently deep into the upper ocean. In the model, heat mainly penetrates to greater depth by slow diffusion processes. In the ocean, however, these transfers are dominated by turbulent mixing. As another consequence of the mixing parameterization, GROB HOPE underestimates mixed-layer depths throughout the year and thus fails to account for Kelvin wave downwelling during El Nino. Thermocline temperatures beneath the mixed layer are about $`2^{}`$ Celsius too warm. Here, the model diffuses too much heat to depths of approximately 500m in the eastern equatorial Pacific which penetrates westward at approximately 250m. This mismatch is the result of unrealistically strong downward diffusion of heat and unrealistically weak upwelling of cold waters.

For Primitive Equation models, vertical transfers pose a greater problem. Nonhydrostatic mixing processes have to be parameterized and such parameterizations are by no means trivial. The mixed-layer model of GROB HOPE is tuned to yield realistic mixing depths at moderate latitudes and compromises for the equatorial mixed layer are accepted. The alternative would be a far more complex and machine-intensive mixed-layer model. Moreover, vertical velocities are determined from mass conservation, independent of the momentum budget. Possible problems and ambiguities are smeared out by diffusion. Hence, models have a tendency to use diffusion where space-time- and phase-space-characteristics of the real ocean are determined by advection and propagation.

Sequential assimilation of subsurface temperatures improves this state estimate significantly. Subsurface temperature data are taken from the TAO/TRITON array which consists of approximately 70 moorings in the tropical Pacific between $`8^{}S`$ and $`8^{}N`$. The buoys record a number of atmospheric parameters, sea surface temperatures and subsurface temperatures at 10 irregularly spaced depths in the upper 500m. Records are transmitted to shore in real-time via the ARGOS satellite system. TAO/TRITON has become one of the most successful ground-based ocean observatories for two major reasons. In the first place, the relatively quiescent tropical waters allow the long-term deployment of buoys. For the more energetic high latitude oceans, the physical lifetime of a similar array would be significantly shorter than time scales of ocean processes of interest. Secondly, the observed variability is readily interpreted in terms of Matsuno’s theory. Since 1985, the TAO/TRITON array has become an integral part of operational services as well as ocean and climate research.

In contrast to the surface boundary condition, nudging of regional TAO/TRITON data is not an option. This array is physically embedded in the global ocean circulation and its data interact physically with their neighborhood in space-time as well as in phase space. In the Fokker-Planck picture of the Kalman Filter this is accounted for by phase space advection and diffusion which are updated at every assimilation step. To this end, model integration is halted at the end of each month and the model temperature field is updated by observed TAO temperatures. After the update, model integration resumes and continues for one month when model temperatures are updated again. Thus, model operation is constrained to the vicinity of the observed state of the ocean.

Fig.5 shows the monthly mean temperature difference Assimilation-Nudge along the equator for December 1997. The data are seen to have three major effects on the estimate: the surface becomes colder, the mixed layer warmer and the thermocline colder. These modifications lead to a significantly higher degree of realism for the estimate. Assimilation ensures that heat supplied at the surface, is uniformly mixed into the upper layer and cold water is upwelled into the thermocline. In response to the data information, the model replaces diffusion dominated dynamics with mixing and advection. However, it is also seen that the assimilation still exhibits some pockets of warm water below 500m although far less than fig.4. Primarily, these pockets are a consequence of the lower boundary condition chosen for vertical transition probabilities: the model is assumed to be true at 500m. Obviously, there is room for improvement.

A different view of these data effects is given in fig.6. The figure shows a time series of monthly mean temperature profiles at a location in the eastern equatorial Pacific for 1997. Simulated (black) and assimilated (red) profiles are compared. The simulation is clearly diffusion dominated, unable to produce a mixed layer and leaking too much heat into the thermocline. The assimilation also fails to produce well-defined mixed layers during the first part of 1997. Before the arrival of the downwelling Kelvin wave in the eastern Pacific, mixed layers here are shallow (typically 25m). Their absence in the assimilation during the first part of the year is a consequence of the poor vertical resolution of GROB HOPE. With the arrival of the Kelvin wave, assimilation produces the characteristic signature of turbulent mixing in the upper ocean with realistic mixed-layer depth. At the same time, the thermocline is cooled by upwelling of colder water and temperature gradients at the mixed-layer base increase.

It is mentioned (corresponding plots not shown) that the model extrapolates data information beyond the temperature field in the TAO/TRITON region. By December 1997, temperatures as far as $`30^{}`$ latitude on both hemispheres are modified by the assimilation. In phase space, model dynamics extrapolate temperature observations onto sealevel, salinity and the velocity field and improve the estimate for the complete state of the ocean. On the other hand, without continued assimilation the model loses memory of the data information from the upper equatorial Pacific after a time period of three to six months.

4.) Summary

Societal needs of climate monitoring call for the installation of global ocean state estimation in the framework of an operational service similar to national weather prediction agencies. Novel ocean observation techniques provide a global data base for such estimates. For a number of parameters, these observations are available almost in real-time and at mesoscale resolution. Assimilators utilize numerical circulation models to dynamically extrapolate measurements in space-time and phase space and, at the same time, constrain model uncertainties.

Given the data volume in ocean observation and modeling, requirements in computing resources are high. For efficiency, it is significant that model and data quality is sufficiently developed to allow nudging the model to essentially unprocessed data at short time constants. This applies particularly to the sea surface where data wealth is largest and determines primarily the surface forcing with a high degree of realism. The forcing problem that plagued ocean model development 20 years ago has been resolved.

Dynamical extrapolation of interior ocean data requires advanced assimilation techniques. Monitoring problems are best addressed by sequential methods such as the Kalman Filter. Typically, observations do not coincide with the corresponding model solution and may even be incompatible with model dynamics. The Kalman Filter finds the initial value for model restart that is both, as close as possible to the data and compatible with model dynamics. To this end, the phase space representation of the filter solves (a large number of) simple, 2+1 dimensional, linear Fokker-Planck equations which represent model dynamics in terms of phase space advection and diffusion. Determination of these parameters by an elementary histogram method circumvents a number of highly complex, but essentially technical issues of the stochastics of nonlinear systems. For numerical analysis, the method proves efficient and reliable.

At this time, models alone are generally unable to simulate the global ocean circulation with satisfactory realism. However, they do capture large-scale features such as gyres, water masses and their seasonal and longer-term variation quite realistically if they are determined by large-scale features of topography and external forcing. Problems typically emerge with the dynamical control of the density field. In the present study, this is demonstrated for equatorial mixing and upwelling and the paths of major ocean currents. On the other hand, the assimilation has shown that models are not antagonistic to (temporary) operation in closer vicinity of the data. At this development stage, model operation in the assimilation mode is capable of delivering practically relevant global ocean state estimates provided a continuous inflow of data is guaranteed.

Some improvement in model performance is readily obtained by fairly simple measures. Higher vertical resolution is oftentimes possible within the framework of given computing resources. Significant increase of the overall horizontal resolution requires generally an upgrade of computer resources. For the performance of HOPE at higher spatio-temporal resolution, see . Moreover, the modeling community continuously develops increasingly appropriate parameterizations of subscale processes and models are updated accordingly.

On the other hand, as long as models are unable to account for basic laws of nature such as the first law of motion, long-term projections with little or no data input after initialization will remain questionable. In this context, it is noted that the “Newtonization” of Richardson’s equations is actually more or less trivial: vertical integration of the Primitive Equations for an incompressible (multilayer) fluid leads again to a Newtonian set of equations of motion. In this framework, vertical variability appears as internal variability of the spatially strictly 2-dimensional fluid. Minor questions arise from the hierarchy problem that results from vertical integration of the nonlinear advection term. Since stratification is represented by the multilayer structure, low-order cut-offs of the advection hierarchy suffice for most purposes of circulation modeling. This shallow water approach to the circulation problem can be formulated with geometric-dynamic integrity and energetic consistency in both, the viscous nonrotating limit and the ideal rotating limit. The theory also admits comprehensive analytical studies of wave-circulation interaction. Currently, physically consistent global circulation models on the basis of shallow water theory are not available.


format: text

text:
# Sensitive Detection of Cold Cesium Molecules by Radiative Feshbach Spectroscopy

## Abstract

We observe the dynamic formation of $`Cs_2`$ molecules near Feshbach resonances in a cold sample of atomic cesium using an external probe beam. This method is $`300`$ times more sensitive than previous atomic collision rate methods, and allows us to detect more than 20 weakly-coupled molecular states, with collisional formation cross sections as small as $`\sigma =3\times 10^{16}`$cm2. We propose a model to describe the atom-molecule coupling, and estimate that more than $`2\times 10^5`$ $`Cs_2`$ molecules coexist in dynamical equilibrium with $`10^8`$ $`Cs`$ atoms in our trap for several seconds.

The newly emerging field of cold molecules provides intriguing possibilities in ultracold gas studies. In the realm of “super-chemistry”, phase-coherent chemical reactions between particles in well-defined internal and external states are dominated by quantum statistics and many-body effects superchem . The coherent nature of the atom-molecule coupling has been recently observed in Bose-Einstein condensates becosci . Quantum manipulation of molecular states has also been proposed for quantum computation qcom . Finally, in precision tests of time-reversal symmetry based on matter-wave interferometry, polar molecules promise a sensitivity several orders of magnitude higher than atoms edmmolecule .

However, the production of cold molecules presents a major experimental challenge. Direct laser cooling of molecules has been proposed using either multiple laser frequencies lasercoolm  or an external cavity cavity ; however, these methods involve exquisite experimental control, and have not yet been implemented. Photoassociation of cold atoms into molecules is limited by the phase space density of the initial atomic sample pa , and consequently only a few thousand molecules have been produced from thermal samples with this method. Two-photon photoassociation from an atomic Bose-Einstein condensate has been limited to similarly small numbers of molecules due to large atom-molecule inelastic collision losses bec2photon . Alternative routes to cold molecules by slowing down a supersonic molecular jet jet , or by buffer gas loading of a magnetic trap doyle , can collect large numbers of molecules, but the resulting temperatures are much higher than in the approaches using cold atoms.

In this work, we employ a radiation-induced loss mechanism to monitor the dynamic formation of cold $`Cs_2`$ molecules near Feshbach resonances. Our method is similar to the detection technique introduced by Heinzen and coworkers heinzen98 ; however, instead of photoassociating ground state atom pairs, which requires knowledge of the excited-state molecular level structure, we use a blue-detuned probe beam to photodissociate them. In a dense sample of optically trapped cesium atoms we observe a class of weakly-coupled $`Cs_2`$ molecular ground states which, in the absence of the probe light, are collisionally metastable at the density of our samples. Since these states are long-lived, a large molecular population can build up even in a thermal sample.

Our method for detecting the formation of ground state molecules, shown in Fig. 1, can be understood as a radiative collision process bur96 ; vul99 . For a typical detuning of $`\mathrm{\Delta }/2\pi +1`$ THz from the free-atom transition, radiative loss occurs when the dipole-dipole interaction $`V_d(r)C_3/r^3`$ between two colliding atoms shifts the excited state into resonance with the laser field at the Condon point $`r_c`$, where $`r_c2`$nm satisfies $`\mathrm{}\mathrm{\Delta }=V_d(r_c)`$ bur96 . The radiative coupling amplitude is therefore proportional to the two-atom wavefunction at the Condon point, $`\mathrm{\Psi }(r_c)`$. On a Feshbach resonance tie9293 , the scattering state $`\mathrm{\Psi }_s(r)`$ couples to a molecular state $`\mathrm{\Psi }_m(r)`$, that has a much larger amplitude at short range and consequently at the Condon point. The coupling amplitude is then enhanced by a factor of $`A|\mathrm{\Psi }_m(r_c)|/|\mathrm{\Psi }_s(r_c)|(r_a/r_m)^{3/2}`$, where $`r_a`$ ($`r_m`$) is the mean interatomic separation in the scattering (molecular) state. Typically, $`r_a1\mu `$m at our atomic density of $`10^{13}`$ cm-3. $`r_m`$ can be estimated from the ground state van der Waals interaction $`V_g=C_6/r^6`$, and the binding energy of the molecular state, which is on the order of the atomic hyperfine splitting $`H_h`$ tie9293 . Given $`C_6=6859`$ a.u. leo00  and $`H_h/h=10GHz`$, we find $`r_m(C_6/H_h)^{1/6}2`$nm. The radiative coupling amplitude of the molecule by the far blue-detuned beam is therefore a factor of $`A10^4`$ stronger than that of two trapped atoms. This indicates that molecules can be detected efficiently without appreciably perturbing the atoms.

Compared to lighter alkalis, cesium atoms display unusually rich cold-collision phenomena, due primarily to the large dipolar interaction ju . The angular dependence of this interaction allows an incident partial wave with angular momentum $`L`$ to couple to molecular states with orbital angular momentum $`N=L\pm 2n\mathrm{}`$, in $`nth`$ order tie9293 . Consequently, many molecular states become available to scattering atom pairs when this interaction is strong. In the following, the above process is referred to as an $`LN`$ resonance, e.g. an $`SG`$ resonance occurs when an s-wave scattering state couples to a molecular state with an orbital angular momentum $`N=4\mathrm{}`$.

In this work, we apply radiative Feshbach spectroscopy to two types of cold collisions: When both scattering atoms are in the lowest hyperfine ground state $`|F=3,m_F=3`$, denoted in the following by $`|3,3+|3,3`$, the scattering state $`\mathrm{\Psi }_s`$ can couple to molecular states $`\mathrm{\Psi }_m`$ with even orbital angular momentum as the relative energy of $`\mathrm{\Psi }_s`$ and $`\mathrm{\Psi }_m`$ is tuned by means of a magnetic field. We find multiple $`SD`$, $`SG`$ Feshbach resonances, and one $`SD`$ shape resonance for samples polarized in $`|3,3`$ fromeite . For $`|3,3+|3,2`$ collisions, on the other hand, both even and odd partial waves are allowed, and we observe an abundance of $`SD`$, $`SG`$, $`PP`$ and $`PF`$ Feshbach resonances fromeite .

Our experimental setup has been described previously in Ref. vul98 . About $`3\times 10^8`$ atoms are loaded into a far-detuned dipole trap formed by a YAG laser at $`1064`$nm using two pulses of Raman-sideband cooling vul98 ; ker00 . We then further optically pump so that up to $`98\%`$ of the remaining atoms are in $`|3,3`$, and by adjusting the optical pumping frequency prepare mixed samples as desired with up to $`10\%`$ in $`|3,2`$. For all of the experiments described here, less than $`1\%`$ ($`0.25\%`$) of the atoms are found to occupy $`|3,1`$ (all other states), as has been determined by microwave spectroscopy chi . Within $`15`$ms the atoms thermalize to a temperature between 3$`\mu `$K and 6$`\mu `$K, and density near $`10^{13}cm^3`$ by means of elastic collisions; the atom number is measured by fluorescence detection and the density is then calculated from the measured temperatures and trap vibration frequencies. We adiabatically ramp the magnetic field from $`0.1`$G to $`2`$G in $`1`$ms, and then to any desired field value up to $`230`$G in another $`1`$ms. No depolarization of the sample is observed during this procedure. A Titanium-Sapphire laser provides the far-detuned probe beam with an $`e^2`$ beam waist of $`2.0`$mm ($`0.6`$mm) in the vertical (horizonal) direction, which uniformly illuminates the sample whose vertical (horizontal) size is $`580\mu `$m ($`60\mu `$m) vul99 . The typical mean intensity and wavelength, optimized for radiative detection sensitivity, are $`20`$W/cm2 and $`847`$nm, respectively. Loss coefficients are extracted from the evolution of the atomic density, assuming a Gaussian distribution, according to $`\dot{n}=n/TGn^2(4/3)^{3/2}Kn^3`$, where $`n`$ is the mean atomic density, $`T=3`$ s is the one-body lifetime, and $`G(K)`$ is the two (three)-body decay coefficient.

Fig. 2 shows the radiative loss rate from $`0`$G to $`230`$G, measured with a resolution of $`100`$mG. All seven observed strong resonances are due to $`|3,3+|3,3`$ collisions. The lower curve in the inset shows a higher-resolution scan performed on a sample with a larger impurity in $`|3,2`$. Seven additional smaller resonances are evident, which originate from $`|3,3+|3,2`$ collisions. The upper curve in the inset is the corresponding measurement of elastic collision-induced evaporative loss performed under similar conditions chi .

It is worth noting that only the strongest resonance at $`48.0`$ G could be detected with the collision rate measurements in Ref. chi . There are two reasons for this: first, there is no appreciable trap loss due to inelastic collisions near these resonances. In particular, for $`|3,3+|3,3`$ collisions, where binary inelastic collisions are endothermic and therefore forbidden, the resonant three-body loss is observed to be below our detection resolution of $`K_3=3\times 10^{27}`$ cm6/s even on the strongest $`SG`$ wave resonance at $`53`$ G. Second, these resonances are sufficiently narrow in energy that the thermally averaged change in elastic collision rate is smaller than our detection limit of $`3\%`$. Clearly the radiative Feshbach resonance detection offers higher sensitivity and accuracy for locating molecular states.

The inset to Fig. 2 shows two complementary aspects of Feshbach resonances: When the field shifts a molecular state into resonance with the scattering continuum, an increase in radiative loss indicates the appearance of a molecular population; at the same time a change in evaporation loss rate indicates a resonant alteration of the effective elastic collision cross section. The difference between these two effects is exemplified by the fact that the radiative resonance at $`47.97`$ G occurs at a slightly lower field value than does the minimum in the evaporation loss at $`48.02`$ G. While the radiative resonance occurs when the molecular state is tuned to the scattering continuum (Feshbach resonance), the minimum in the elastic cross-section occurs when the background scattering length is cancelled by the resonant contribution from the bound state (Ramsauer-Townsend effect). The finite width of the resonance causes these two points to be distinct. For stronger (wider) resonances this is more dramatic: the wide minimum in the elastic cross section observed at $`17.1`$G in Ref.chi  is in fact induced by the Feshbach resonance at a negative field of $`8.25`$G (or equivalently, positive field for the $`|3,3`$ state) fromeite .

We model the radiative process using a chemical rate-equation approach that treats atoms and molecules as two different species timmermans . In a thermal sample of de Broglie wavelength $`\lambda _{dB}`$ with low phase space density $`\varphi =n\lambda _{dB}^31`$, the evolution of the mean atomic (molecular) density $`n`$ ($`m`$) near a Feshbach resonance can be described classically: Based on detailed balance between one molecular state and $`1/\varphi `$ scattering states in the relative coordinate, $`m=\frac{1}{2}\varphi n`$ molecules per unit volume can coexist with $`n`$ atoms in thermal equilibrium in the absence of collision loss nphi . The formation and destruction rates of molecules are then described by $`\alpha \varphi n/2`$ and $`\alpha m`$, where $`\alpha `$ characterizes the strength of the Feshbach resonance. Given a leading order collisional loss $`\gamma nm`$ due to atom-molecule collisions, and radiative atom loss induced by the probe beam $`L=2\delta mI`$, where $`\gamma `$ and $`\delta `$ are rate constants, and $`I`$ is the probe beam intensity, we obtain the following rate equations:

$`\dot{m}`$ $`=`$ $`\alpha \varphi n/2\alpha m\gamma nmL/2`$ (1)
$`L`$ $`=`$ $`2\delta Im.`$ (2)
The stationary solution for the molecular density $`m`$ and the fractional atom loss rate $`L/n`$ are then given by

$`m`$ $`=`$ $`{\displaystyle \frac{(\varphi /2)n}{1+(\gamma /\alpha )n+(\delta /\alpha )I}}`$ (3)
$`L/n`$ $`=`$ $`{\displaystyle \frac{\delta I\varphi }{1+(\gamma /\alpha )n+(\delta /\alpha )I}}.`$ (4)
To test this model, we measure the radiative atom loss as a function of probe beam intensity $`I`$ for different atomic densities, as shown in Fig. 4 for the strongest $`SG`$ resonance at $`53.50`$ G. Here the mean atomic density is switched between $`n_1=1.2\times 10^{13}`$cm-3 and $`n_2=0.6\times 10^{13}`$cm-3 by changing the trap loading parameters. The phase space densities are $`\varphi _11/100`$ and $`\varphi _21/200`$, respectively. The loss constant $`L/n`$ is derived from the measured radiation exposure time $`t`$ necessary to deplete $`10\%`$ of the sample, namely, $`L/n=0.1/t`$. The polarization purity of $`95\%`$ in $`|3,3`$ and the temperature of $`5.5\mu `$K do not vary systematically by more than $`5\%`$ between the two settings. We have verified that no additional loss is due to field ramping or radiative processes induced by the very far detuned dipole trap at $`1064`$nm.

At low probe beam intensity $`I<5`$W/cm2, the fractional loss rate is linear in $`I`$. This suggests that the molecular and atomic populations reach a dynamic equilibrium where the probe beam only weakly depletes the molecules, without appreciably changing their density. At high probe beam intensity $`I>10`$W/cm2, the fractional loss rate is proportional to the sample density and is independent of the probe beam intensity. This implies a quick depletion of molecular population by the probe beam, leading to a loss rate that is therefore equal to the formation rate of the molecules. As indicated by the solid lines in Fig. 4, Eq. 4 agrees well with the experimental results, and a simultaneous fit to both sets of data with only two free parameters yields $`\alpha =364(20)`$s-1 and $`\delta /\gamma =3.3(6)\times 10^{12}`$(Wcm)-1

From the above measurements, we can first deduce the sensitivity of radiative Feshbach spectroscopy for detecting molecule formation. Given $`\alpha `$, the molecule formation cross section at $`53`$ G is $`\sigma _\alpha =\alpha \varphi /vn=8\times 10^{14}`$cm2, where $`v`$ is the mean velocity of the atoms. While this resonance induces a maximum loss rate of $`4`$s-1, the loss rate of the weakest observed resonance at $`7.8`$G is approximately $`0.04`$s-1, with a signal-to-noise ratio of $`3`$. We estimate that the sensitivity of our detection, expressed in terms of a molecule formation cross section, is then $`\sigma _\alpha /3003\times 10^{16}`$cm2. Compared to the nominal elastic cross section for atoms in $`|3,3`$ of $`1.5\times 10^{10}`$cm2, this sensitivity corresponds to the creation of less than two molecules every $`10^6`$ collisions.

Furthermore, we can calculate the stationary number of molecules in the gas from the measured rate constants. While the latter only reveal the formation constant $`\alpha `$, we can safely put an upper limit on the collisional loss $`\gamma `$ from the unitarity limit of the s-wave cross section for inelastic collisions: $`\sigma <\pi /k^2`$. We then estimate $`\gamma =\sigma _\gamma v_\gamma <(\pi /k^2)(\mathrm{}k/\mu )`$, where $`\mathrm{}^2k^2/2\mu \frac{3}{2}k_BT`$ is the collision energy, $`\mu =\frac{2}{3}M`$ the reduced mass, and $`M`$ the mass of the Cs atom. Given Eq. 3 and $`\alpha `$, we obtain $`\gamma <2.8\times 10^{11}cm^3`$s-1 and a ratio of molecular to atomic density $`m/n>2.6(5)\times 10^3`$ at the atomic phase space density of $`\varphi =1/100`$ and low radiative excitation. Equivalently, more than $`2.6(5)\times 10^5`$ molecules coexist with $`10^8`$ atoms in our trap.

In conclusion, we have used radiative Feshbach spectroscopy to detect more than twenty new $`Cs_2`$ molecular states with high sensitivity and resolution. Since no strong collision loss accompanies these resonances, as many as $`3\times 10^5`$ molecules accumulate and coexist dynamically with $`10^8`$ atoms in our dipole trap. In the future, we hope to explore the regime where the atom-molecule coupling strength exceeds the thermal energy of the gas, and its coherent nature may become evident even in our thermal sample.

This work was supported in parts by grants from AFOSR and the NSF.


text:
# Origin of the anomalies: The modified Heisenberg equation

## Abstract

The origin of the anomalies is analyzed. It is shown that they are due to the fact that the generators of the symmetry do not leave invariant the domain of definition of the Hamiltonian and then a term, normally forgotten in the Heisenberg equation, gives an extra contribution responsible for the non conservation of the charges. This explanation is equivalent to that of the Fujikawa in the path integral formalism. Finally, this formalism is applied to the conformal symmetry breaking in two-dimensional quantum mechanics.

The use of the symmetries of a system is one of the most fruitful techniques in physics, specially for quantum systems. Among the consequences of symmetry in quantum physics are: selection rules, relations between matrix elements of observables, degeneracies in energy and, specially, the existence of conservation laws which are guarantied by the Noether’s theorem or by its equivalent, the Ehrenfest equation. In particular, the evolution of the expectation values of an operator $`B`$ is given by the Heisenberg equation

$`{\displaystyle \frac{d}{dt}}<\mathrm{\Psi }(t)|B\mathrm{\Psi }(t)>=`$ $`<\mathrm{\Psi }(t)|{\displaystyle \frac{B}{t}}\mathrm{\Psi }(t)>+{\displaystyle \frac{i}{\mathrm{}}}<\mathrm{\Psi }(t)|[H,B]\mathrm{\Psi }(t)>,`$ (1)
that says that for any group of symmetry whose elements (or generators in the associated Lie algebra) commute with the Hamiltonian $`H`$, if such operators do not depend explicitly on time $`t`$, their expectation value on any physical state must be constant with $`t`$.

Even for the case when the symmetry is not exact (it is explicitly broken), the Heisenberg equation (1) gives us how do evolve with time the expectation values of the corresponding generators. This has been largely used, for example, in particle physics where the flavor symmetry is explicitly broken by mass terms, or in nuclear physics where isospin symmetry is broken by the Coulomb interaction between protons and by the up-down quark mass difference.

However, there are some cases where although the symmetry is exact at the classical level, it is not preserved in the corresponding quantum theory. This is the anomalous symmetry breaking that was first discovered in quantum field theory when studding certain Feynman diagrams, and further in the analysis of $`\pi ^02\gamma `$ decay ABJ , also in the Schwinger model Sch . In the path integral formalism, the existence of anomalies can be viewed as a consequence of that, in this case, even if the classical Lagrangian is invariant under the symmetry, the measure is not FUY , that gives rise to extra surface terms which originate the anomaly.

In the Hamiltonian formalism, the anomaly can be understood as a consequence from the fact that the Heisenberg equation (1) is exact only if the domain of definition of the Hamiltonian is invariant by the operator $`B`$; in other cases, it appears an extra term which is the responsible for the anomaly JGE . To be more precise, let $`H`$ be the quantum Hamiltonian which is self-adjoint when defined on a domain $`D_H`$, then for any physical state in the Hilbert space $``$ and any operator $`B`$ we have

$`{\displaystyle \frac{d}{dt}}<\mathrm{\Psi }(t)|B\mathrm{\Psi }(t)>`$ $`=`$ $`<\mathrm{\Psi }(t)|{\displaystyle \frac{B}{t}}\mathrm{\Psi }(t)>+{\displaystyle \frac{i}{\mathrm{}}}\left(<H\mathrm{\Psi }(t)|B\mathrm{\Psi }(t)><\mathrm{\Psi }(t)|BH\mathrm{\Psi }(t)>\right)`$
which can be written as

$`{\displaystyle \frac{d}{dt}}<\mathrm{\Psi }(t)|B\mathrm{\Psi }(t)>=<\mathrm{\Psi }(t)|{\displaystyle \frac{B}{t}}\mathrm{\Psi }(t)>+{\displaystyle \frac{i}{\mathrm{}}}<\mathrm{\Psi }(t)|[H,B]\mathrm{\Psi }(t)>+𝒜,`$ (2)
where the anomaly $`𝒜`$ is defined as

$`𝒜={\displaystyle \frac{i}{\mathrm{}}}<\mathrm{\Psi }(t)|\left(H^+H\right)B\mathrm{\Psi }(t)>,`$ (3)
and in (2) the commutator $`[H,B]`$ is defined in the whole Hilbert space. Equivalently, in the Heisenberg picture, we obtain that for any operator $`B`$ their derivative with respect to the time obeys a generalized Heisenberg equation:

$`{\displaystyle \frac{dB}{dt}}={\displaystyle \frac{B}{t}}+{\displaystyle \frac{i}{\mathrm{}}}[H,B]+{\displaystyle \frac{i}{\mathrm{}}}(H^+H)B.`$ (4)
Comparing equations (1) and (2) we see that the Heisenberg equation (1) is exact only whenever $`B`$ keeps invariant the domain of definition $`D_H`$ of the Hamiltonian, because if $`B\mathrm{\Psi }_hD_H`$ for $`\mathrm{\Psi }_hD_H`$, then the extra term $`<\mathrm{\Psi }_h(t)|\left(H^+H\right)B\mathrm{\Psi }_h(t)>=0`$ as long as $`H^+=H`$ when acting on states of $`D_H`$. But when $`B`$ does not keep $`D_H`$ invariant, that is $`B\mathrm{\Psi }_hD_H`$, the extra term will give a surface contribution responsible for the anomaly. In general, it is said that in the presence of an anomaly, the commutator of the Hamiltonian with the corresponding charges has two contributions, the regular one and an extra part originated by the anomaly, $`[H,B]_{total}=[H,B]_{reg}+[H,B]_{anom}`$. We see that the so called regular part is nothing but the commutator of the extension of the operators to the whole Hilbert space and the anomalous term is just $`(H^+H)B`$.

It can be proved that the above description of the anomalies is equivalent to that of the path integral FUY . For example, in quantum mechanics, the Feynman propagator is

$`K(z,t;y,0)={\displaystyle _{\stackrel{x\left(0\right)=y}{x(t)=z}}}[d\mu (x)]\mathrm{exp}(iS(x,\dot{x})/\mathrm{})={\displaystyle _n\phi _n^{}(y)\phi _n(z)e^{(iE_nt/\mathrm{})}},`$ (5)
where $`\phi _n`$ are the eigenvectors of the Hamiltonian and $`_n`$ means sum over the discrete and integral over the continuum spectrum. In this sense, different self-adjoint extensions $`H^{(\lambda )}`$ of the quantum Hamiltonian associated to the same classical Lagrangian, give rise to different sets of orthonormal eigenvalues $`\phi _n^{(\lambda )}(x)`$, depending on the self-adjoint extension defined on $`D_H^{(\lambda )}`$, and each of them is characterized by a different measure in the path integral version, so if a particular domain of definition $`D_H^{(\lambda )}`$ is not invariant under the operator $`B`$, the same is true for the associated measure. The proof for quantum field theories is equivalent.

It should be noted that the existence of an anomaly is independent of the need or not of a regularization process for the theory, as can be seen in some quantum mechanical systems as that of a charged particle moving on a two-torus and coupled to an electromagnetic field, see NSM  and also AAE .

Recently, there has been a renewed interest in anomalies in conformal quantum mechanics: the 3-dimensional $`1/r^2`$ potential which is relevant as an example of an anomaly in Molecular Physics CGC , and the two dimensional $`\delta `$ interaction RJ . In what follows, we shall analyze the later on the light of the generalized Heisenberg equation (2). The problem is that of a free particle in two-dimensional quantum mechanics with a $`\delta ^2(r)`$ interaction.

$`H={\displaystyle \frac{P^2}{2m}}+\lambda {\displaystyle \frac{1}{r}}\delta (r).`$ (6)
It can be seen that considering the extension of $`H`$ to the whole Hilbert space $`=L^2(R_+,rdr)L^2(S_1,d\phi )`$, the theory is scale invariant, that is, the dilation operator $`D=tHG`$ (where $`G=(1/4)(\mathrm{𝐫𝐩}+\mathrm{𝐩𝐫})`$), the conformal generator $`K=t^2H+2tD+(1/2)r^2`$ and $`H`$ close on commutation:

$`{\displaystyle \frac{i}{\mathrm{}}}[K,D]`$ $`=`$ $`K,`$ (7)
$`{\displaystyle \frac{i}{\mathrm{}}}[H,K]`$ $`=`$ $`2D,`$ (8)
$`{\displaystyle \frac{i}{\mathrm{}}}[D,H]`$ $`=`$ $`H,`$ (9)
showing that the invariance algebra is $`SO(2,1)`$. The above equations together with the classical Heisenberg equation (1) mean that $`\frac{d}{dt}<K>=\frac{d}{dt}<D>=0`$, and that there can not be any normalizable bound state. However, in order to properly define the quantum theory, we must first define the domain of definition $`D_H`$ of the Hamiltonian in such a way that $`\overline{D}_H=`$, and that when acting on $`D_H`$ we have $`H^+=H`$. In order to do that, we start by removing the origin to avoid the singularities of the Hamiltonian, then working on $`\dot{𝐑}^\mathrm{𝟐}=𝐑^\mathrm{𝟐}/\{\mathrm{𝟎},\mathrm{𝟎}\}`$ AGH  and in polar coordinates, so the Hamiltonian reads:

$`H={\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{^2}{r^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{^2}{\phi ^2}}\right).`$ (10)
The first question is to define the domain of definition of the operator $`d^2/d\phi ^2`$ which has deficiency indices $`d_+=d_{}=2`$, so there are infinitely many self-adjoint extensions associated to different physical situations . Here we shall start with the particular one with periodic boundary conditions, so defining $`d^2/d\phi ^2`$ on

$`d_0=\left\{f(\phi )L^2(S_1,d\phi )\right|f(0)=f(2\pi );f^{}(0)=f^{}(2\pi )\},`$ (11)
in that case, acting on $`d_0`$ the eigenfunctions of $`d^2/d\phi ^2`$ are $`\xi _n(\phi )=(2\pi )^{1/2}e^{in\phi }`$ with $`nZ`$, and then $`D_H`$ can be written as

$`D_H={\displaystyle \underset{nZ}{}}\left(D_n(R_+,rdr){\displaystyle \xi _n(\phi )}\right),`$ (12)
where $`D_n(R_+,rdr)`$ must be chosen in such a way that the radial part of the Hamiltonian

$`H_r={\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{d^2}{dr^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}{\displaystyle \frac{n^2}{r^2}}\right),`$ (13)
is self-adjoint. For $`n0`$ the deficiency indices of $`H_r`$ are (0,0) so it is essentially self-adjoint acting on functions which vanishes at $`r=0`$ and infinity. But for $`n=0`$, the deficiency indices are $`(1,1)`$ and there are infinitely many self-adjoint extensions characterized by a parameter $`\beta `$:

$`D_{n=0}^\beta =\left\{f(r)L^2(R_+,rdr)\right|\underset{r0}{lim}\left({\displaystyle \frac{f(r)}{\mathrm{log}(\alpha _0r)}}\right)=\beta \underset{r0}{lim}\left[f(r)\underset{r^{}0}{lim}\left({\displaystyle \frac{f(r^{})}{\mathrm{log}(\alpha _0r^{})}}\right)\mathrm{log}(\alpha _0r)\right]\},`$
(14)
where $`\alpha _0^2=\left(2m\mathrm{\Lambda }_0/\mathrm{}^2\right)`$, and $`\mathrm{\Lambda }_0`$ is the dimensional constant we must introduce in order to make the deficiency index equations $`H_r\mathrm{\Psi }_\pm (r)=\pm i\mathrm{\Lambda }_0\mathrm{\Psi }_\pm (r)`$ dimensionally consistent. The case $`\beta =0`$ is the Friedrich’s extension and corresponds to the situation $`\lambda =0`$ in (6), whereas $`\beta 0`$ accounts for the case $`\lambda 0`$.

Equation (14) means that if $`f(r)D_{n=0}^\beta `$ then for $`r0;f(r)a(\mathrm{log}(\alpha r)+b)+\theta (r)`$ with

$`{\displaystyle \frac{1}{\beta }}=b+\mathrm{log}({\displaystyle \frac{\alpha }{\alpha _0}}).`$ (15)
Now, it can be seen that for $`n0`$, we have $`GD_nD_n`$ so $`G`$ leaves invariant $`D_{n0}`$, but for $`n=0,\beta 0`$, and $`f(r)D_0^\beta `$, we obtain

$`Gf(r)D_0^\beta ^{};\beta ^{}={\displaystyle \frac{\beta }{\beta +1}},`$ (16)
so $`D_0^\beta `$ is not invariant by the action of $`G`$, and consequently it is not invariant by $`D`$ and $`K`$. Hence the symmetry will be anomalously broken. The most relevant manifestation of this anomalous symmetry breakdown is the existence of a normalizable bound state

$`\mathrm{\Psi }_0(r,\phi )={\displaystyle \frac{\alpha }{\pi ^{1/2}}}K_0(\alpha r),`$ (17)
with energy

$`E_0={\displaystyle \frac{\mathrm{}^2}{2m}}\alpha ^2,`$ (18)
where if $`\mathrm{\Psi }_0D_0^\beta `$ then $`\alpha `$ and $`\beta `$ are related by $`1/\beta =\mathrm{log}(\alpha /2\alpha _0)+\gamma `$, ($`\gamma `$ is the Euler’s constant). In this case for the dilation operator $`D`$, the left side of equation (2) evaluates to $`E_0`$ whereas for the right side we have that $`<\mathrm{\Psi }_0|\frac{D}{t}\mathrm{\Psi }_0>=\frac{i}{\mathrm{}}<\mathrm{\Psi }_0|[H,D]\mathrm{\Psi }_0>`$, and the equation (2) reduces to

$`{\displaystyle \frac{d}{dt}}<\mathrm{\Psi }_0|D\mathrm{\Psi }_0>=𝒜={\displaystyle \frac{i}{\mathrm{}}}<(H^+H)D\mathrm{\Psi }_0>.`$ (19)
Integrating by parts, we finally obtain

$`𝒜`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2\alpha ^2}{2m}}\left\{\left[r{\displaystyle \frac{d}{dr}}K_0(\alpha r)(r{\displaystyle \frac{d}{dr}}+2)K_0(\alpha r)\right]_0^{\mathrm{}}\left[rK_0(\alpha r){\displaystyle \frac{d}{dr}}(r{\displaystyle \frac{d}{dr}}+2)K_0(\alpha r)\right]_0^{\mathrm{}}\right\}`$ (20)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2\alpha ^2}{2m}}.`$
Then we see that it is precisely the fact that $`D`$ does not keep $`D_0^\beta `$ invariant which originates the extra contribution which accounts for the value of the left side.

Finally, we can consider what happens with other self-adjoint extensions of the Hamiltonian. For example, defining the operator $`d^2/d\phi ^2`$ with vanishing boundary conditions on $`d_v=\{f(\phi )L^2(S_1,d\phi )|f(0)=f(2\pi )=0\}`$, everything, respective to the anomaly, is equivalent to the case of periodic boundary conditions; but if we define $`d^2/d\phi ^2`$ on $`d_\theta =\{f(\phi )L^2(S_1,d\phi )|f(0)=e^{i2\pi \theta }f(2\pi );f^{}(0)=e^{i2\pi \theta }f^{}(2\pi );\theta [0,1)\}`$ then the domain of definition of $`H`$ can be written as $`D_H=_{nZ}(D_{n,\theta }(R_+,rdr)(2\pi )^{(1/2)}e^{i(n+\theta )\phi })`$, and for $`\theta 0`$ the deficiency indices of the radial part of the Hamiltonian

$`H_r^\theta ={\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{d^2}{dr^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}{\displaystyle \frac{(n+\theta )^2}{r^2}}\right),`$ (21)
are $`(0,0)`$ for $`n0,1`$, and $`(1,1)`$ for $`n=0,1`$. In this last case, defining the adequate domains $`D_{0,\theta }^\beta `$ and $`D_{1,\theta }^\beta ^{}`$ it is easy to see that both subspaces are not invariant by $`G`$, which results in the fact that there are two normalizable bound states.

In conclusion, it has been shown that the origin of the anomalous symmetry breakdown is that the generators of the symmetry do not leave invariant the domain of definition of the Hamiltonian and then, although the formal commutator of those generators with $`H`$ vanishes, the charges are not conserved due to the extra surface term that appears in the exact Heisenberg equation (2). For the case of the conformal symmetry breaking in the $`\delta ^2(r)`$ potential, the anomaly has been calculated exactly and the Eq. (2) verified. Similar results for the $`1/r^2`$ potential in three dimensional quantum mechanics will be discussed elsewhere.

Acknowledgments

I would like to thank M. Asorey, for useful discussions. This work has been supported the CICYT grants BFM2000-1057 and FPA2000-1252.


text:
# Transport anisotropy in biaxially strained La2/3Ca1/3MnO3 thin films

(received June 04, 2002)

## Abstract

Due to the complex interplay of magnetic, structural, electronic, and orbital degrees of freedom, biaxial strain is known to play an essential role in the doped manganites. For coherently strained La2/3Ca1/3MnO3 thin films grown on SrTiO3 substrates, we measured the magnetotransport properties both parallel and perpendicular to the substrate and found an anomaly of the electrical transport properties. Whereas metallic behavior is found within the plane of biaxial strain, for transport perpendicular to this plane an insulating behavior and non-linear current-voltage characteristics (IVCs) are observed. The most natural explanation of this anisotropy is a strain induced transition from an orbitally disordered ferromagnetic state to an orbitally ordered state associated with antiferromagnetic stacking of ferromagnetic manganese oxide planes.

to appear in Physical Review B (Brief report)

It is well known that the physics of the doped perovskite manganites is determined by a complex interplay of structural, magnetic, electronic, and orbital degrees of freedom. While the classical double exchange model can qualitatively explain the transition from a paramagnetic insulating to a ferromagnetic metallic state Zener:51 , for a more complete understanding of the physics of the manganites electron-lattice coupling has to be included Millis:95 . Recently, Millis et al. have pointed out that uniform compression, as realized by hydrostatic pressure, increases the electron hopping amplitude favoring a ferromagnetic metallic state Millis:98 . In contrast, biaxial strain, as realized in epitaxial thin films grown on substrates with significant lattice mismatch, enhances the Jahn-Teller distortions favoring an insulating state due to the tendency of the electrons to become localized Millis:98 . Fang et al. have calculated the phase diagram of the almost tetragonal doped manganite La1-xSrxMnO3 as a function of biaxial strain by studying the instabilities of the ferromagnetic state Fang:00 . Their results are in ageement with experiments on biaxially strained La1-xSrxMnO3 thin films Konishi:99 . In their work, it has been shown that the orbitally disordered ferromagnetic state (F) is unstable against orbital ordered states with layer (A) and chain (C) type antiferromagnetism. In turn, it is expected that these different magnetic states are associated with different magnetotransport behavior via the double exchange mechanism. The further investigation of the validity of the strain phase diagram, the corresponding magnetotransport properties, and the extendability to doped manganites with strong tilt of MnO6 octahedra Vigliante:2001  and phenomena as charge ordering as for example La1-xCaxMnO3 is of great interest to gain more insight into the physics of these materials and its dependence on lattice distortions. For the purpose of this study, it is important to verify the coherency of the strained state of the doped manganite, in order to be able to determine properly the intrinsic properties within the biaxial strain phase diagram, and to exclude effects of nonuniform strain distribution or relaxation effects Biswas:00 ; Biswas:01 .

In this Letter we present a careful study of the structural, electronic, and magnetic properties of coherently strained La2/3Ca1/3MnO3 (LCMO) thin films and LCMO-La2/3Ba1/3MnO3 (LBMO) heterostructures on SrTiO3 substrates. In addition to previous studies, transport properties have been measured both parallel and perpendicular to the plane of biaxial strain. The key result of our study is that biaxial strain results in highly anisotropic transport properties: Whereas insulating behavior and non-linear IVCs are observed perpendicular to the biaxially strained plane (parallel to the $`c`$ axis), the in-plane transport is metallic below the Curie temperature $`T_C`$. The saturation magnetization of biaxially strained LCMO is strongly reduced compared to the bulk material or the less strained LBMO films. It is shown that this behavior is not due to interface effects between different layers Bibes:01 , but is an intrinsic property of the biaxially strained LCMO arising most likely from a strain induced orbital ordering. Another important result is that in strained LCMO a low-resistance state can be induced by applying either a magnetic field or a high current density. A similar behavior is for example observed in Pr0.7Ca0.3MnO3 and Nd0.5Sr0.5MnO3 and has been attributed to a magnetic field or current induced (local) melting of a charge resp. orbitally ordered ground state Asamitsu:97 ; Fiebig:98 ; Guha:00a .

We have grown LCMO films and LBMO-LCMO-LBMO heterostructures on SrTiO3 ($`a3.905`$ Å) substrates using pulsed laser deposition Gross:2000a ; Klein:2001a . The lattice mismatch between the SrTiO3 substrate and LCMO ($`a_{\mathrm{bulk}}3.864`$ Å in pseudocubic notation) is -1.2% resulting in in-plane tensile strain, whereas the mismatch between SrTiO3 and LBMO ($`a_{\mathrm{bulk}}3.910`$ Å) is only about 0.13% resulting in very small compressive strain. Note that in the LBMO-LCMO-LBMO heterostructures the LBMO layers provide low resistance contacts to the (ultra)thin LCMO films allowing for a homogeneous current feed for transport perpendicular to the film. Furthermore, effects of a possible surface “dead layer” Sun:99  are avoided in the multilayer structure. Layer-by-layer growth of the films was confirmed by a high pressure reflection high energy electron diffraction (RHEED) system showing clear growth oscillations Rijnders:97 ; Klein:00 .

As has already been stressed it is important to prove the coherency of the strained state. The coherent film thickness was determined from Laue oscillations in $`\theta 2\theta `$ x-ray scans. It was found that both LCMO and LBMO films grow coherently strained on SrTiO3 substrates up to a thickness of at least 60 nm, consistent with literature Bibes:01 ; Izumi:98 . The out-of-plane ($`c`$ axis) and in-plane ($`a`$ axis) lattice parameters of the LCMO thin films grown on SrTiO3 have been determined from the (002) and (103) reflections. The in-plane film lattice parameters coincide throughout the whole layer with the substrate lattice parameters. That is, the in-plane lattice parameter of LCMO is enlarged, while the out-of-plane lattice constant is reduced, leading to a tetragonal lattice distortion with $`c/a0.985`$. The tetragonal distortion can be viewed as a Jahn-Teller like distortion resulting in an increased Jahn-Teller splitting of the Mn $`e_g`$ levels and, in turn, in a tendency of the electrons to become localized.

In order to measure the electrical transport properties perpendicular to the film plane, mesa structures (see inset of Fig. 1) were patterned using optical lithography and Ar ion beam etching. The mesas have typical area of several $`\mu `$m2. Fig. 1 shows the resistance vs temperature, $`R(T)`$, curves for LBMO-LCMO-LBMO mesa structures with different thickness $`d_{\mathrm{LCMO}}`$ of the LCMO layer. For comparison, a sample without LCMO layer ($`d_{\mathrm{LCMO}}=0`$) is shown. It is evident that the resistance increases with increasing $`d_{\mathrm{LCMO}}`$. We estimate the involved resistivities $`\rho _c`$ to about several $`\mathrm{\Omega }`$m at 10 K. Note that due to the non-linear IVCs, it is difficult to obtain a meaningful $`\rho _c`$ or thickness dependence $`\rho _c(d)`$. In the inset we show the sheet resistance as a function of $`d_{\mathrm{LCMO}}`$. The strong increase of resistance with increasing layer thickness can partially be due to tunneling through the barrier. The important point is that the resistivity does not show any tendency to saturation meaning the resistance is intrinsic to the barrier and not to the interface. Furthermore, the metallic $`R(T)`$ behavior below $`T_C`$ that is observed for the sample with $`d_{\mathrm{LCMO}}=0`$ turns into an semiconducting or insulating $`R(T)`$ behavior with increasing $`d_{\mathrm{LCMO}}`$. The low-temperature plateau (below 60 K) of the resistivity corresponds to a plateau in the temperature dependence of the magnetization (see upper panel of Fig. 3) as expected for double exchange materials. This thickness dependence clearly shows that the insulating $`R(T)`$ behavior is not due to the patterning process, due to a contact resistance between the gold contact layer and the manganite film, nor due to interface effects. One can conclude that the insulating behavior for transport along the $`c`$ axis observed at low temperatures is an intrinsic property of the coherently strained LCMO thin films. We note that in the work of Bibes et al. interfaces between magnetic LCMO and non-magnetic SrTiO3 are shown to favor phase segregation Bibes:01 . Inhomogeneous transport properties were found in La0.7Ca0.3MnO3/SrTiO3 heterostructures and interpreted as arising from magnetic interface disorder Jo:99 . In our case, the interfaces are between two differently doped manganites. While it remains to be investigated whether interdiffusion is present at such interfaces, the effects observed here are related to the magnetic order of the whole thin film layer.

In Fig. 2 we show the IVCs of several mesa structures with different values of $`d_{\mathrm{LCMO}}`$ measured in $`c`$ axis direction, i. e. with current perpendicular to the biaxially strained plane. While for low enough current densities $`VI`$, all mesas show highly non-linear IVCs for higher current densities following a $`VI^n`$ dependence with $`n0.20.3`$. In contrast, the IVCs measured with the current in-plane are ohmic. For comparison, mesa structures have been patterned into single LBMO films ($`d_{\mathrm{LCMO}}=0`$) that are well lattice matched to the SrTiO3 substrate. For these samples, ohmic IVCs have been obtained both for the current applied in- and out-of-plane. From this and the fact that the measured voltage clearly increases with increasing $`d_{\mathrm{LCMO}}`$ (see discussion before), preparation or interface effects can be excluded as the origin of the nonlinear IVCs. Hence, we can conclude that for current perpendicular to the plane of biaxial strain the insulating behavior of the coherently strained LCMO films is associated with non-linear current transport. We note that the non-linearity becomes smaller with increasing $`T`$ and vanishes at the $`T_C`$ of the LCMO thin films ranging between 100 and 150 K. $`T_C`$ was determined from magnetization measurements of the bi/trilayer films before mesa patterning. This strongly suggests that the electronic anisotropy is coupled to the magnetic properties of the LCMO films.

In agreement with previous experiments we found a strain dependence of the saturation magnetization $`M_S`$ of the LCMO films. Zandbergen et al. reported $`M_S2.5\mu _B`$/Mn atom at 5 K for a 6 nm thick, coherently strained LCMO film on SrTiO3 Zandbergen:99 . Consistently, we have measured a value of $`M_S2.2\mu _B`$/Mn atom at 5 K for a 57.5 nm thick LCMO film (see inset of Fig. 2). In contrast, almost strain free LBMO films on SrTiO3 show the expected saturation magnetization of about $`3.67\mu _B`$/Mn atom. These findings are in agreement with the phase diagram predicted by Fang et al. Fang:00  where for tensile strain an instability to an A-type antiferromagnetic state is predicted. Another interesting observation is that the size of the hysteresis loop in $`M(H)`$ curves also depends on strain. While for almost strain free LBMO films on SrTiO3 a coercive field of $`\mu _0H_c10`$ mT is observed at 5 K, a much larger value of $`\mu _0H_c70`$ mT is measured for the strained LCMO films.

Fig. 3 gives an overview on the field, current and temperature dependence of the resistance and magnetization in strained LCMO films by showing $`R(T)`$ curves recorded at different applied fields and currents. Applying high magnetic fields results in a strong suppression of the resistance at all $`T`$. Also, due to the non-linearity of the IVCs the measured resistance is reduced when the applied current is increased below about 150 K. We note that for $`T150`$ K the resistance is dominated by the LCMO layer, only around 250 K it is dominated by the LBMO layer. Fig. 3 clearly shows that with increasing field, the non-linearity becomes weaker and also the onset temperature of the non-linear behavior is shifted from about 100 K at 0 T to 50 K at 8 T. For comparison, the inset of Fig. 3 shows the in-plane $`\rho (T)`$ curves of a LCMO film ($`d_{\mathrm{LCMO}}=57.5`$ nm). Clearly, a metallic $`\rho (T)`$ behavior is observed below the peak temperature. The same film shows an insulating $`R(T)`$ behavior for current perpendicular to plane.

To further study the effect of the applied current, we have measured IVCs for successive current sweeps with increasing amplitude as shown in Fig. 4. After zero field cooling, the current is first increased to 10 $`\mu `$A (curve a) corresponding to a current density of about 100 A/cm2. On decreasing the current again (curve b), a lower voltage is measured at the same current values, i.e. the applied current of 10 $`\mu `$A has switched the sample to a state with lower resistance. Applying a high magnetic field (8 T) has the same effect as applying a high current (1 mA). After applying a field of 8 T, the measured IVC are stable and the resistance is independent on the applied current. Thus, both a high magnetic field and a high current density induce a state with reduced resistivity.

We now address the possibility of a phase separated state in fully strained LCMO. We note that in the case of inhomogeneously strained or relaxed films (e.g. due to island growth), it is plausible to assume a phase separated state with ferromagnetically and antiferromagnetically ordered clusters, as has been discussed recently for LCMO thin films grown on LaAlO3 Biswas:00 ; Biswas:01 . A similar phenomenon has been observed for bulk Pr0.7Ca0.3MnO3 samples Fiebig:98  as well as for LCMO bulk and thin film samples Faeth:99 . However, phase separation cannot explain the transport anisotropy present in our samples, but would be expected to lead to direction independent behavior.

All our experimental observation can be naturally described by the assumption of strain induced orbital order as predicted by Fang et al. Fang:00 . For tensile strain ($`c/a<1`$, as present in our samples), a transition from the conventional double exchange mediated, orbital disordered F state to the orbital ordered A state, which is composed mainly by $`d_{x^2y^2}`$ states, is expected. Whereas in the F state the spins are aligned parallel in adjacent planes, in the A state anti-parallel orientation of the ferromagnetically ordered planes is present. That is, the gradual transition from the F to the A state is accompanied by a strong reduction of saturation magnetization in agreement with the strongly reduced value obtained in our experiments. Furthermore, the A-type antiferromagnetic state can be metallic only within the ferromagnetically ordered plane, but is insulating in perpendicular direction. This again is in agreement with our experimental observation. That is our data can be interpreted by a strain induced orbital ordering effect at fixed doping. Evidently, a sufficiently high current density allows to switch between the competing states with high current density or high field favoring the F state which can be interpreted as (local) melting of the orbital order by introduction of highly spin-polarized carriers.

In summary, we have investigated coherently strained LCMO films. The biaxial strain was found to induce anisotropic transport properties at low temperatures with metallic and insulating behavior for current in- and out-of-plane. It has been shown that this behavior is not due to inhomogeneous interface effects or phase separation. We suggest strain induced orbital ordering as the origin of the observed behavior in agreement with theoretical predictions. We also have shown that by applying a high current density, a low-resistance state can be induced. This effect may be of interest with respect to magnetoelectronic devices.

This work was supported by the Deutsche Forschungsgemeinschaft. One of us, J. K., acknowledges support by the Graduiertenkolleg 549.