text:
# Quantum Circa-Rhythms

## I Introduction

Perhaps the greatest achievement of a pre-quantum phenomenology in chemistry was the Mendeleev periodic table of elements. The table, in particular, allowed to predict new elements with correctly anticipated properties. But to really understand the structure of the table one had to invent a simply-looking equation

$`i|\dot{\psi }=H|\psi `$ (1)
where $`|\psi `$ is a vector from a Hilbert space and $`H`$ a self-adjoint operator.

The Schrödinger equation (1) with appropriately chosen $`H`$ allowed to understand the origin and types of chemical bonds, that is, the deep structures behind the Mendeleev table. Although the identification of $`H`$ for a concrete molecule was of crucial importance, there were certain features of the abstract equation that allowed to determine general properties of those systems whose dynamics could be governed by (1). For example, a self-adjoint $`H`$ was related to conservative and isolated sytems. The conserved quantity is the average energy $`E=\psi |H|\psi `$. Conservative systems that are non-isolated (i.e. non-trivially interact with the outside world) were described by another simply-looking equation

$`i\dot{\rho }=[H,\rho ]`$ (2)
where $`\rho `$ is a trace-class positive operator acting in the same Hilbert space and $`H`$ is the same self-adjoint operator. Eq. (2) is known as the von Neumann equation. (2) is equivalent to (1) if $`\rho `$ projects on the direction of $`|\psi `$. The invariant $`\psi |H|\psi `$ is replaced by $`E=\mathrm{Tr}(H\rho )`$. If $`|j`$ is a normalized vector, then the diagonal element $`\rho _{jj}=j|\rho |j=\mathrm{Tr}P_j\rho =p_j`$ is a probability associated with the proposition represented by $`P_j=|jj|`$. The probability is time dependent if $`[H,P_j]0`$.

Thinking of the present-day chemistry or biochemistry one is struck by the ubiquity of nonlinear kinetic equations Prigogine ; Goldbeter . Physically, all kinetic equations describe dynamics of probabilities and in this sense are close to (2). Typical kinetic equations describe open systems and the associated probabilities are related with commuting propositions. Nonlinearity occurs whenever there are feedbacks or cooperative phenomena, and details of evolution depend on the number of substrates and phenomenology used to model right-hand-sides of the equations.

Similarly to the catastrophe theory, phenomenological modelling is to some extent an art Thom . Having experience with nonlinear evolutions one is capable of choosing the right-hand-sides of the equations in such a way that certain qualitative elements of the dynamics will have to occur, and then it remains to play with parameters. Not trying to diminish the role of phenomenological modelling we have to admit that phenomenology is a two-edged sword. On the one hand, it links qualitative analysis with quantitative predictions, and often reveals aspects that are initially unexpected. On the other hand, however, shapes of the equations are dictated by our needs. This should be contrasted with (1) and (2) which can be derived from first principles.

A general property of systems described by (1) or (2) is the exact linearity of the evolutions they generate. Linear evolutions practically do not occur in biological and chemical systems, so linearity is in this context a drawback. Sometimes a linear approximation is justifed if some paramters are small, but the generic case is nonlinear. In Ref. I  we have argued that a candidate equation describing a conservative quantum system in a feedback with some environment is the nonlinear generalization of the von Neumann equation

$`i\dot{\rho }=[H,f(\rho )]`$ (3)
where $`f(\rho )`$ is an operator function satisfying the “no-feedback, no-nonlinearity” property, formally meaning that $`f(\rho )=\rho `$ if $`\rho `$ is a projector on some direction $`|\psi `$ in a Hilbert space. Eq. (3) is integrable in the sense of soliton theory for any $`f`$. The time-invariant is now $`E=\mathrm{Tr}\left(Hf(\rho )\right)`$. A soliton technique of solving (3), based on Darboux transformations, was introduced in LC , further developed in UCKL , and recently generalized in CCU . Therefore, as opposed to the kinetic equations that typically have to be solved numerically (for exceptions cf. num1 ; num2 ; num3 ), we can work with exact analytic solutions. This allows us to have access to the time scales and precision numerically unavailable.

The set of solutions of (3) contains all the solutions of form $`\rho =|\psi \psi |`$ where $`|\psi `$ satisfies (1). A new class of “self-switching” solutions occuring only for nonlinear $`f`$ was discovered in LC  and applied in I  in the context of morphogenesis. These solutions are characterized by three regimes involving a single switching between two asymptotically linear evolutions.

Still, looking for the most striking manifestations of biochemical nonlinear evolutions one arrives at phenomena that are rhythmic. Periodic switching, with periods ranging from fraction of a second to years, is encountered at all levels of biological organization, and its successful modelling in terms of nonlinear kinetic equations is one of the great achievements of computational biology (for a review cf. Goldbeter2002 ).

The kinetic equations considered in phenomenological models of biochemical oscillations involve nonlinearities of a mixed linear-Hill type. The dynamical systems are typically finite dimensional, the number of degrees of freedom corresponding to the number of different substrates one takes into account. The use of Hill functions of the form $`\chi _{y,n}(x)=x^n/(y^n+x^n)`$ is motivated by phenomenological arguments and their ability to continuously interpolate between 0 and 1. The latter feature allows for a dynamical mechanism of switching, a property whose origin can be traced back to bifurcations of nonlinear maps and catastrophe theory.

Of particular interest are theoretical studies of circa-rhythms. Let us recall that circa-rhythms are “classes of rhythms that are capable of free-running in constant conditions with a period approximating that of the environmental cycle to which they are normally synchronized” Moore . The examples of the rhythms are the circadian (24 hours), the circalunar (28 hours), or the circannual (365.25 days) ones. What makes circa-rhythms interesting from a formal point of view is the interplay between the phases where there exists an external forcing (light–darkness periods, say), and the rhythmicity that sustains even after the external driving is switched off (as in experiments in constant darkness). Formally, the two cases may correspond to nonlinearities involving (or not) an explicit time dependence of some parameters.

To give an example, the model of circadian rhythms in Neurospora Gonze2000  involves four variables and two types of Hill functions. The rhythmicity occurs in the model even without an external driving, but the case with explicit time dependence of coefficients is treated as well. Similar constructions are given in Gonze2002a ; Gonze2002b .

In the present paper we address the problem of circa-rhythmicity from the perspective of Eq. (3). The main idea is to consider a composite Hamiltonian system that, as a whole, is conservative. There are two sets of degrees of freedom corresponding to its two subsystems. When we look at certain averages associated with one of the subsystems we find that their time evolution can be determined by an effective density matrix whose dynamics is given again by an equation of the form (3), but now with explicitly time-dependent coefficients in the nonlinear function $`f(\rho )`$. Therefore, the interaction between the two subsystems effectively turns one of them into a pacemaker that drives the other subsystem and generates its evolution with characteristic bursting patterns exibiting “day” and “night” phases. As opposed to the phenomenological models found in the literature we do not “force” switching by putting Hill functions by hand. We take a simple quadratic nonlinearity satisfying the “no feedback, no nonlinearity” condition. The motivation for this concrete choice comes from its simplicity and experience showing that the nonlinear behavior found here is generic and analogous to the one found for more complicated $`f`$’s (cf. the example of $`f(\rho )=\rho ^q2\rho ^{q1}`$ with arbitrary real $`q`$ analyzed in UCKL ). The solutions we obtain are also interesting in themseleves as first examples of solutions of (3) exhibiting more than one self-swiching.

The paper is organized as follows. We begin with an example of a nonlinear von Neumann equation (3) for a composite system and time independent nonlinearity. We then find its particular solution that is reducible to a new von Neumann equation for a subsystem but with explicitly time-dependent nonlinearity. We then generalize the example to the class of time-dependent von Neumann equations and show how to construct their solutions. Next we explicitly solve an example and illustrate its rhythmic properties.

## II Preliminaries: Dynamics in a subsystem

As a simple example consider

$`i\dot{\varrho }`$ $`=`$ $`[\omega (\underset{J}{\underset{}{b_1^{}b_1b_2^{}b_2}})+\left(1+X\right)a^{}a,\varrho ][Xa^{}a,f(\varrho )]`$ (4)
where $`[b_k,b_k^{}]=1=[a,a^{}]`$ and

$`X=b_1+b_1^{}+b_2+b_2^{}`$ (5)
For $`\varrho ^2=\varrho `$ we find

$`i\dot{\varrho }`$ $`=`$ $`[\omega J+a^{}a,\varrho ]`$ (6)
describing three independent degrees of freedom: harmonic oscillation with unit frequency combined with rotation in a plane with frequency $`\omega `$. For $`\varrho ^2\varrho `$ rotation and oscillation get nonlinearly coupled.

Let us now eliminate the rotation by switching to a rotating reference frame

$`\rho =e^{i\omega Jt}\varrho e^{i\omega Jt}`$ (7)
where the density matrix satisfies

$`i\dot{\rho }`$ $`=`$ $`[\left(1+X(t)\right)a^{}a,\rho ][X(t)a^{}a,f(\rho )].`$ (8)
Denoting $`Y=i(b_1+b_2^{}b_1^{}b_2)`$ we find

$`X(t)`$ $`=`$ $`e^{i\omega Jt}Xe^{i\omega Jt}=X\mathrm{cos}\omega t+Y\mathrm{sin}\omega t`$
$`Y(t)`$ $`=`$ $`e^{i\omega Jt}Ye^{i\omega Jt}=Y\mathrm{cos}\omega tX\mathrm{sin}\omega t`$
and, for any $`t`$ and $`t^{}`$,

$`[X(t),X(t^{})]=[Y(t),Y(t^{})]=[X(t),Y(t^{})]=0.`$ (9)
Since $`X`$, $`Y`$, and $`a^{}a`$ commute we can introduce their joint eigenvectors

$`X|x,y,n`$ $`=`$ $`x|x,y,n`$ (10)
$`Y|x,y,n`$ $`=`$ $`y|x,y,n`$ (11)
$`a^{}a|x,y,n`$ $`=`$ $`n|x,y,n`$ (12)
$`X(t)|x,y,n`$ $`=`$ $`x(t)|x,y,n`$ (13)
$`Y(t)|x,y,n`$ $`=`$ $`y(t)|x,y,n`$ (14)
$`x(t)`$ $`=`$ $`x\mathrm{cos}\omega t+y\mathrm{sin}\omega t`$ (15)
$`y(t)`$ $`=`$ $`y\mathrm{cos}\omega tx\mathrm{sin}\omega t.`$ (16)
Now take any time-independent normalized vector

$$|\psi =𝑑x𝑑y\psi (x,y)|x,y$$

and make the Ansatz $`\rho (t)=|\psi \psi |w(t)`$ where $`w`$ is a density matrix acting only on the oscillator degress of freedom. Denoting $`g(\rho )=\rho f(\rho )=|\psi \psi |g(w)`$ we rewrite (8) as

$`|\psi \psi |i\dot{w}=|\psi \psi |[a^{}a,w]`$ (17)
$`+X(t)|\psi \psi |a^{}ag(w)|\psi \psi |X(t)g(w)a^{}a`$
Taking matrix elements of both sides of (17) between arbitrary $`x,y|`$ and $`|x^{},y^{}`$ we obtain

$`\psi (x,y)\overline{\psi (x^{},y^{})}i\dot{w}=\psi (x,y)\overline{\psi (x^{},y^{})}[a^{}a,w]`$ (18)
$`+\psi (x,y)\overline{\psi (x^{},y^{})}\left(x(t)a^{}ag(w)x^{}(t)g(w)a^{}a\right).`$
The Ansatz is internally consistent only if $`\psi (x,y)=\sqrt{\delta (xx_0)\delta (yy_0)}`$. The diagonal $`x=x^{}=x_0`$, $`y=y^{}=y_0`$ leads to

$`i\dot{w}`$ $`=`$ $`[a^{}a,w]x_0(t)[a^{}a,g(w)],`$ (19)
$`x_0(t)=x_0\mathrm{cos}\omega t+y_0\mathrm{sin}\omega t`$, which is an equation of the form (3) but with $`f`$ replaced by

$`\stackrel{~}{f}(w)=wx_0(t)g(w)=\left(1x_0(t)\right)w+x_0(t)f(w).`$ (20)
If $`f`$ satisfies the “no feedback, no nonlinearity” condition $`f(\rho )=\rho `$ for $`\rho ^2=\rho `$, the same holds for the effective $`\stackrel{~}{f}`$.

## III Quantum pacemaker and feedback

The pacemaker is an oscillatory system whose state can be modified by an external entraining agent, a zeitgeber Moore . Our main interest in this paper is in the free-running oscillation of the pacemaker $`x_0(t)`$, and the circa-rhythms it induces in the absence of the zeitgeber.

The next level is the coupling between the pacemaker, which defines the clock “mechanism”, and the observable-level circa-rhythms which define the “hands” of the clock. We assume a nonlinear feedback between the hands and the mechanism, but the dynamics is nondissipative in the sense that the energy of the hands averaged over a single cycle of the oscillation is constant.

Motivated by the analysis of the previous section we concentrate on the following class of equations

$`i\dot{w}(t)=f_1(t)[H,w(t)]+f_2(t)[H,f\left(w(t)\right)]`$ (21)
which define the state of the hands of the clock. The “no-feedback, no-nonlinearity” condition reads $`w(t)=f_1(t)w(t)+f_2(t)f\left(w(t)\right)`$ whenever $`w(t)^2=w(t)`$. Now assume we know a solution $`w_0(t)`$ of (3). Then

$`w(t)`$ $`=`$ $`e^{iH_0^tf_1(x)𝑑x}w_0\left(_0^tf_2(x)𝑑x\right)e^{iH_0^tf_1(x)𝑑x}`$ (22)
is a solution of (21) as can be verified by a direct calculation. The whole problem of solving (21) reduces to finding a solution of (3), which can be performed by soliton techniques LC ; UCKL .

Let us consider a simple but generic example where $`f(w)=(1s)w+sw^2`$ and

$`f_1(t)=1+ϵ\mathrm{cos}\omega t,f_2(t)=ϵ\mathrm{cos}\omega t`$ (23)
The parameter $`s`$ allows us to compare situations where the driven dynamics is linear ($`s=0`$) and purely nonlinear ($`s=1`$), and for the two cases investigate the role of the $`ϵ`$s. Varying $`s`$ we can also investigate stability properties of the rhythms under fluctuations of the feedback.

## IV Hands of the clock

The hands of the clock are described by the Hamiltonian $`H=a^{}a=_{n=0}^{\mathrm{}}n|nn|`$ of a harmonic oscillator type. The frequency of the oscillator is equal to unity, meaning that this is a reference frequency for all the other frequencies found in the system. Modelling the hands by a harmonic oscillator is quite natural for obvious reasons. The quantum oscillator has in addition the appealing property of being delocalized. The oscillations occur at the level of probabilities in position space and what we call “hands” is a kind of a center-of-mass quantity. This type of oscillation is what one expects in a biological system, since classical oscillators such as pendula or springs are too “mechanical” to be of any relevance. Intuitively, in our model, the hands move at time $`t`$ to the region of space where the concentration of probability is the biggest. The probability density may be regarded as a measure of state of a certain spatially extended object. In the absence of a feedback between the hands and the pacemaker, the hands oscillate harmonically with their own internal frequency. We shall see later that the nonlinear coupling may practically suppress the internal oscillations of the hands in certain intervals of time. What will remain are the sudden bursts occuring, roughly, with the period of the pacemaker.

Quantum harmonic oscillator is an infinite dimensional dynamical system and, hence, a solution of von Neumann equations may be characterized by an arbitrary number of parameters determining the initial state of the hands. There exists a simple trick allowing to construct infinitely-dimensional solutions on the basis of a single finite-dimensional one. The trick exploits equal spacing of eigenvalues of $`H`$. With our choice of units the eigenvalues are given simply by natural numbers. Let us divide them into sets containing $`N`$ elements: $`\{0,1,\mathrm{},N1\}`$, $`\{N,N+1,\mathrm{},N+N1\}`$, and so on. Each such subset corresponds to a block in the Hamiltonian, and each block can be represented by a $`N\times N`$ diagonal matrix of the form

$`H_k=k\mathrm{𝟏}+\mathrm{diag}(0,1,\mathrm{},N1)=k\mathrm{𝟏}+H_0.`$ (24)
As a consequence, in each block we have to solve the same matrix equation since a restriction $`w_k`$ of $`w`$ to the $`k`$-th subspace satisfies

$`i\dot{w}_k=[H_k,f(w_k)]=[H_0,f(w_k)].`$ (25)
The job can be reduced to finding a sufficiently general solution of a $`N\times N`$ problem. In each subspace we can take a different initial condition and a different normalization of trace. The whole infinite-dimensional solution will take the form of a direct sum

$`w(t)=_{k=0}^{\mathrm{}}p_kw_k(t,p_k),`$ (26)
$`_{k=0}^{\mathrm{}}p_k=1`$. The $`k`$-th part depends on $`p_k`$ in a complicated way since the function $`f(w)`$ is not 1-homogeneous, i.e. $`f(\lambda w)\lambda f(w)`$. The inhomogeneity implies that change of normalization simultaneously rescales time; the normalization of probability implies that a change of $`p_k`$ in a $`k`$-th subspace influences all the other subspaces by making their dynamics faster or slower. In this sense the solution, in spite of its simplicity, is not a simple direct sum of independent evolutions.

In order to illustrate the possible effects we can use the solutions derived in I  for the simplest nontrivial case involving self-switching, i.e. for $`N=3`$ and quadratic nonlinearity. Of course, since $`N=3`$, the Hamiltonian $`H`$ and the solution $`w`$ possess matrix elements describing transitions $`12`$, $`23`$, $`13`$, between three different basis vectors. Although the observation is in itself trivial, we want to stress here the formal analogy to the so-called three-variable models of bursting oscillations in enzyme reactions with autocatalytic regulation DecrolyGoldbeter1982 ; DecrolyGoldbeter1987 .

We select a subspace spanned by three subsequent vectors $`|k`$, $`|k+1`$, $`|k+2`$. The family of interest is parametrized by $`\alpha 𝐑`$ controlling the “moment” and type of switching between bursts. The parameter naturally occurs at the level of the Darboux transformation, where it characterizes an initial condition for the solution of the Lax pair. The density matrices $`w_k(t)=_{m,n=0}^2w_{mn}|k+mk+n|`$ are completely characterized by the $`k`$-independent matrix of time-dependent coefficients $`w_{mn}`$. The reader may check by a straightforward substitution Math  that the matrix

$`\left(\begin{array}{ccc}w_{00}& w_{01}& w_{02}\\ w_{10}& w_{11}& w_{12}\\ w_{20}& w_{21}& w_{22}\end{array}\right)={\displaystyle \frac{1}{15+\sqrt{5}}}\left(\begin{array}{ccc}5& \xi (t)& \zeta (t)\\ \overline{\xi }(t)& 5+\sqrt{5}& \xi (t)\\ \overline{\zeta }(t)& \overline{\xi }(t)& 5\end{array}\right)`$ (33)
with

$`\xi (t)`$ $`=`$ $`{\displaystyle \frac{\left(2+3i\sqrt{5}i\right)\sqrt{3+\sqrt{5}}\alpha }{\sqrt{3}\left(e^{\gamma t}+\alpha ^2e^{\gamma t}\right)}}e^{i\omega _0t},`$
$`\zeta (t)`$ $`=`$ $`{\displaystyle \frac{9e^{2\gamma t}+\left(1+4\sqrt{5}i\right)\alpha ^2}{3\left(e^{2\gamma t}+\alpha ^2\right)}}e^{2i\omega _0t}`$
is indeed a normalized ($`\mathrm{Tr}w=1`$) solution of the von Neumann equation

$`i\dot{w}=[H,(1s)w+sw^2],`$ (34)
$`H=\mathrm{diag}(0,1,2)`$. The parameters are $`\omega _0=1\frac{5+\sqrt{5}}{15+\sqrt{5}}s`$, $`\gamma =\frac{2}{15+\sqrt{5}}s`$.

Now let us rescale the trace. We do it in three steps. The modified density matrix

$`w_1(t)=e^{i(1s)Ht}w(t)e^{i(1s)Ht}`$ (35)
is a solution of $`i\dot{w}_1(t)=[sH,w_1(t)^2].`$ Therefore

$`w_2(t)=\mathrm{\Lambda }e^{i(1s)\mathrm{\Lambda }Ht}w(\mathrm{\Lambda }t)e^{i(1s)\mathrm{\Lambda }Ht}`$ (36)
is also a solution of $`i\dot{w}_2(t)=[H,sw_2(t)^2]`$ and

$`w_3(t)=e^{i(1s)Ht}w_2(t)e^{i(1s)Ht}`$ (37)
is a solution of

$`i\dot{w}_3(t)=[H,(1s)w_3(t)+sw_3(t)^2].`$ (38)
Performing these operations on our explicit solution we find

$`w_3(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{15+\sqrt{5}}}\left(\begin{array}{ccc}5& \xi _3(t)& \zeta _3(t)\\ \overline{\xi }_3(t)& 5+\sqrt{5}& \xi _3(t)\\ \overline{\zeta }_3(t)& \overline{\xi }_3(t)& 5\end{array}\right)`$ (42)
with

$`\xi _3(t)`$ $`=`$ $`{\displaystyle \frac{\left(2+3i\sqrt{5}i\right)\sqrt{3+\sqrt{5}}\alpha }{\sqrt{3}\left(e^{\gamma \mathrm{\Lambda }t}+\alpha ^2e^{\gamma \mathrm{\Lambda }t}\right)}}e^{i\left(\omega _0\mathrm{\Lambda }+(1s)(1\mathrm{\Lambda })\right)t},`$
$`\zeta _3(t)`$ $`=`$ $`{\displaystyle \frac{9e^{2\gamma \mathrm{\Lambda }t}+\left(1+4\sqrt{5}i\right)\alpha ^2}{3\left(e^{2\gamma \mathrm{\Lambda }t}+\alpha ^2\right)}}e^{2i\left(\omega _0\mathrm{\Lambda }+(1s)(1\mathrm{\Lambda })\right)t},`$
Now $`_0^tf_1(x)𝑑x=t+\frac{ϵ}{\omega }\mathrm{sin}\omega t`$, $`_0^tf_2(x)𝑑x=\frac{ϵ}{\omega }\mathrm{sin}\omega t`$. We finally obtain the solution

$`w(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{15+\sqrt{5}}}\left(\begin{array}{ccc}5& \xi (t)& \zeta (t)\\ \overline{\xi }(t)& 5+\sqrt{5}& \xi (t)\\ \overline{\zeta }(t)& \overline{\xi }(t)& 5\end{array}\right)`$ (46)
with

$`\xi (t)`$ $`=`$ $`{\displaystyle \frac{\left(2+3i\sqrt{5}i\right)\sqrt{3+\sqrt{5}}\alpha }{\sqrt{3}\left(e^{\gamma \mathrm{\Lambda }\frac{ϵ}{\omega }\mathrm{sin}\omega t}+\alpha ^2e^{\gamma \mathrm{\Lambda }\frac{ϵ}{\omega }\mathrm{sin}\omega t}\right)}}`$
$`\times e^{i\left(ϵ\omega _0\mathrm{\Lambda }+ϵ(1s)(1\mathrm{\Lambda })ϵ\right)\frac{\mathrm{sin}\omega t}{\omega }}e^{it}`$
$`\zeta (t)`$ $`=`$ $`{\displaystyle \frac{9e^{2\gamma \mathrm{\Lambda }\frac{ϵ}{\omega }\mathrm{sin}\omega t}+\left(1+4\sqrt{5}i\right)\alpha ^2}{3\left(e^{2\gamma \mathrm{\Lambda }\frac{ϵ}{\omega }\mathrm{sin}\omega t}+\alpha ^2\right)}}`$
$`\times e^{2i\left(ϵ\omega _0\mathrm{\Lambda }+ϵ(1s)(1\mathrm{\Lambda })ϵ\right)\frac{\mathrm{sin}\omega t}{\omega }}e^{2it}`$
of

$`i\dot{w}(t)`$ $`=`$ $`\left(1+ϵ\mathrm{cos}\omega t\right)[H,w(t)]`$ (47)
$`ϵ\mathrm{cos}\omega t[H,(1s)w(t)+sw(t)^2]`$
Since $`\mathrm{Tr}Hw(t)=\mathrm{Tr}Hw(0)`$, $`\mathrm{Tr}Hw(t)^2=\mathrm{Tr}Hw(0)^2`$, the internal energy of the system, averaged over one period $`T`$ of the pacemaker oscillation, is

$`E={\displaystyle \frac{1}{T}}{\displaystyle _t^{t+T}}𝑑t^{}\mathrm{Tr}Hf\left(w(t^{})\right)=\mathrm{Tr}Hw(0),`$ (48)
and does not depend on $`t`$. In this sense the subsystem is conservative.

If one does not integrate over the pacemaker period, one finds that the internal energy of the hands harmonically oscillates with the pacemaker frequency $`\omega `$. What is characteristic, however, the hands do not oscillate harmonically but behave as if they were accumulating energy during the phases of quiescence in order to suddenly release it in violent bursts. Moreover, the bursts are different during the two halves of the pacemeaker cycle, and thus resemble the day-night differences one finds in real organisms.

All these properties are illustrated on explicit examples. We plot the dynamics of the hands $`\overline{q}(t)=\mathrm{Tr}\widehat{q}w(t)=_{\mathrm{}}^{\mathrm{}}𝑑qqq|w(t)|q`$ as functions of time for different parameters characterizing the nonlinearity, and for different initial conditions.

## V Discussion

The ability of plants and mammals to measure approximately 24-hour intervals of time is no doubt of an endogenous origin. Our current understanding of circadian rhythms leads to a notion of a pacemaker, an oscillating timing system whose role is to coordinate various complicated and different in shape rhythms of a biological system.

Splitting of a “clock” into a “pacemaker” and “hands” occurs in our example quite naturally. The entire “clock” is a conservative system so the rhythmic behavior we find is of endogenous origin. However, the system cannot be treated as isolated from external environment since the solution $`\rho =Pw`$ is not a projector and, hence, involves nontrivial external correlations (cf. the analysis of this point given in I ).

To make our analysis simple we have played with a concrete family of exact solutions of (3). The family is characterized by several parameters related to the choice of the equation ($`\omega `$, $`s`$) and the initial condition for the dynamics ($`x_0=ϵ`$, $`\alpha `$). The solutions we have derived reveal an unknown aspect of nonlinear von Neumann equations — the possibilty of evolution involving an infinite number of switching events.

In spite of its simplicity, the quadratic nonlinearity generates circa-rhythmic behavior of surprising complexity and variety. Much more complicated types of behavior would have been found if one had used the freedom of higher dimensionality offered by (26). Moreover, the solution of the factorized form $`Pw`$, with $`P`$ a projector, is very special and occurs only in cases the operators $`X(t)`$ and $`X(t^{})`$ commute for all $`t`$, $`t^{}`$. Had we replaced the Abelian group of rotations in a plane (with generator $`J`$) by something more complicated, we would not have been able to use the Ansatz with factorization.

To end these remarks let us make it clear that we do not claim to have found a better description of nonlinear oscillations occuring in biochemical systems than the phenomenological approach of Goldbeter and his Brussels school. What we do claim, however, is that nonlinear von Neumann equations of the form (3), in addition to their appealing general properties, imply a vast range of phenomena with analogies in biophysical systems. One should not be surprised if (3) will one day play a fundamental role in life sciences.

###### Acknowledgements.

Our work was supported by the KBN grant No. 5 P03B 040 20 (MC) and the Flemish Fund for Scientific Research (FWO project G.0335.02).


format: text

text:
# The origin and tidal evolution of cuspy triaxial haloes

## 1 INTRODUCTION

A generic prediction of the currently favoured cold dark matter (CDM) cosmological model of hierarchical structure formation is that dark matter (DM) haloes of galaxies and clusters are flattened triaxial systems (e.g., Barnes & Efstathiou, 1987; Frenk et al., 1988; Dubinski & Carlberg, 1991; Warren et al., 1992; Cole & Lacey, 1996; Thomas et al., 1998; Jing & Suto, 2002). Observationally, inferring the intrinsic shape of DM haloes of galaxies and clusters is a difficult task, but nonetheless several potentially powerful probes exist that may allow us to to distinguish between spherical and flattened DM haloes. These include the dynamical modelling of collisionless tracers such as tidal streams orbiting the Milky Way (Johnston et al., 1999; Ibata et al., 2001; Mayer et al., 2002; Majewski et al., 2004), the distribution and kinematics of gas in spiral galaxies (Kuijken & Tremaine, 1994; Franx, van Gorkom & de Zeeuw, 1994; Schoenmakers, Franx & de Zeeuw, 1997; Merrifield, 2002), and polar ring galaxies (Schweizer, Whitmore & Rubin, 1983; Sackett & Sparke, 1990; Sparke, 2002; Iodice et al., 2003) which provide shape constraints perpendicular to the disc plane, gravitational lensing applications (Kochanek, 1995; Bartelmann, Steinmetz & Weiss, 1995; Koopmans et al., 1998; Oguri, Lee & Suto, 2003) and the flattening of the extended X-ray isophotes in elliptical galaxies (Buote & Canizares, 1994, 1996, 1998; Buote et al., 2002).

Numerical experiments of isolated equilibrium models are very useful for studying the dynamical evolution of gravitating systems in a controlled way. Several techniques exist for constructing $`N`$-body realizations of spherical haloes and multi-component galaxies (Hernquist, 1993; Holley-Bockelmann et al., 2003; Kazantzidis, Magorrian & Moore, 2004a, hereafter KMM), but it is far more complicated to build triaxial equilibria. Even though the Jeans theorem guarantees that the distribution function depends only on the isolating integrals of motion, explicit expressions for the latter other than the energy per unit mass, $`E`$, are rarely known in the case of triaxial potentials. As a result, the triaxial models that have been constructed so far have been either limited to a few special analytical cases (e.g., Stäckel potentials or rotating $`f(E_J)`$ models with $`E_J`$ being the Jacobi constant) or entirely based on numerical techniques (Schwarzschild, 1979, 1993; Merritt & Fridman, 1996; Poon & Merritt, 2001; Terzić, 2003).

Recently, Holley-Bockelmann et al. (2001) used a technique of adiabatically applying a drag to the velocities of the particles along each principal axis in order to create cuspy triaxial systems starting with a spherical Hernquist (Hernquist, 1990) model. The original prescription of Hernquist (1993) was generalized by Boily, Kroupa & Peñarrubia-Garrido (2001) for accomodating composite, axisymmetric models of galaxies and extended by Tinker & Ryden (2002) for investigating the effect of rotating, triaxial halos on disk galaxies. These techniques are useful and have been used to understand the response of systems to live triaxial potentials. However, they are somewhat restricted to modest values of the flattening and it is difficult to incorporate a significant amount of rotational flattening.

Triaxiality can arise in a number of ways. The most common way of constructing numerical models of triaxial galaxies is via the merger of other objects. Examples include binary mergers of spherical haloes (e.g., White, 1978; Fulton & Barnes, 2001) and disk galaxies (e.g., Gerhard, 1981; Barnes, 1992; Barnes & Hernquist, 1996; Naab & Burkert, 2003; Kazantzidis et al., 2004b), as well as multiple mergers of systems (e.g., Weil & Hernquist, 1996; Dubinski, 1998). The structure of the final remnant of two component models depends sensitively on both the orbital geometry (Naab & Burkert, 2003) as well as on the inclination and internal properties of the disks (Kazantzidis et al., 2004b).

A general technique for constructing cuspy axisymmetric and triaxial $`N`$-body systems would be advantageous for many purposes including studying the dynamical friction in flattened and rotating systems, the mass loss and tidal stripping from triaxial substructure haloes, the effects of baryonic accretion and disc formation on triaxial halo shape and anisotropy, the effects of triaxial shape on the properties and stability of discs, the interaction between central black holes and cusps and the inflow of gas in triaxial systems.

In this paper we explore the generation of triaxial structures that is based on the way that haloes in hierarchical models of structure formation obtain their shapes. Inspection of a cosmological simulation reveals that triaxiality arises via mergers that take place with various amounts of angular momentum which is generated from the large-scale tidal field. Mergers between two haloes that occur with little angular momentum (radial mergers) produce prolate systems, whilst high angular momentum mergers produce oblate systems. Most CDM haloes form via a sequence of mergers with varying amounts of angular momentum such that haloes with arbitrary triaxiality may be formed and the triaxiality will vary with radius (Moore et al., 2001; Vitvitska et al., 2002).

The outline of this paper is as follows. In Section 2, we describe our technique for constructing equilibrium cuspy $`N`$-body haloes with various degrees of flattening. There we discuss in detail the numerical experiments we performed and present our results for the internal structure of the resulting models. In Section 3, we investigate the tidal evolution of triaxial substructure haloes within a static host potential and study the shape of subhaloes within a cosmological CDM simulation. Finally, we summarise our main results in Section 4.

## 2 CONSTRUCTING CUSPY TRIAXIAL DARK MATTER HALOES

We begin with a cuspy spherical DM halo which follows the Navarro, Frenk & White (1996, hereafter NFW) density profile

$$\rho (r)=\frac{\rho _\mathrm{s}}{(r/r_\mathrm{s})(1+r/r_\mathrm{s})^2}\text{(}rr_{\mathrm{vir}}\text{)}.$$
(1)
Here the characteristic inner density $`\rho _\mathrm{s}`$ and scale radius $`r_\mathrm{s}`$, are sensitive to the epoch of halo formation and tightly correlated with the halo virial parameters, via the concentration, $`c`$, and the virial overdensity $`\mathrm{\Delta }_{\mathrm{vir}}`$. Since the NFW density profile corresponds to a cumulative mass distribution that diverges as $`r\mathrm{}`$, we impliment an exponential cut-off for $`r>r_{\mathrm{rvir}}`$ which sets in at the virial radius and turns off the profile on a scale $`r_{\mathrm{decay}}`$ which is a free parameter and controls the sharpness of the transition

$$\rho (r)=\frac{\rho _\mathrm{s}}{c(1+c)^2}\left(\frac{r}{r_{\mathrm{vir}}}\right)^ϵ\mathrm{exp}\left[\frac{rr_{\mathrm{vir}}}{r_{\mathrm{decay}}}\right](r>r_{\mathrm{vir}}),$$
(2)
where $`cr_{\mathrm{vir}}/r_\mathrm{s}`$ is the concentration parameter. Finally, in order to ensure a smooth transition between (1) and (2) at $`r_{\mathrm{vir}}`$, we require the logarithmic slope there to be continuous. This implies

$$ϵ=\frac{1+3c}{1+c}+\frac{r_{\mathrm{vir}}}{r_{\mathrm{decay}}}.$$
(3)
Once the density profile $`\rho (r)`$ of our initial model has been specified, we construct an $`N`$-body realization using the technique discussed in detail in KMM. To summarize, the particle positions and velocities are initialized by sampling the exact phase-space distribution function using Eddington’s inversion formula (Eddington, 1916) for the given density profile $`\rho (r)`$ and assuming spherical symmetry. Note that in this paper we will consider only initial models with isotropic velocity dispersion tensors and therefore their distribution function will only depend on the binding energy $`E`$, $`f=f(E)`$.

For the following experiments we begin with an initial galaxy-sized halo with virial mass $`M_{\mathrm{vir}}=7\times 10^{11}h^1\mathrm{M}_{}`$ and concentration of $`c=10`$ in accordance with theoretical predictions for objects at this mass scale and for the adopted cosmology (e.g., Bullock et al., 2001). Note that for all simulations presented in this paper including the cosmological ones described in the next section, we consider as our framework the concordance flat $`\mathrm{\Lambda }`$CDM model with present-day matter and vacuum densities $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, respectively, dimensionless Hubble constant $`h=0.7`$, present-day fluctuation amplitude $`\sigma _8=0.9`$ and index of the primordial power spectrum $`n=1.0`$. For this choice of concentration, the virial and scale radius of the initial halo model are $`r_{\mathrm{vir}}=181h^1\mathrm{kpc}`$ and $`r_\mathrm{s}=18.1h^1\mathrm{kpc}`$, respectively. The initial model is constructed with $`N=5\times 10^5`$ particles and we use a gravitational softening of $`ϵ=0.3h^1\mathrm{kpc}`$.

All the simulations presented in this paper were carried out using PKDGRAV, a multi-stepping, parallel $`N`$-body code (Stadel, 2001) which uses a spline kernel softening and multi-stepping based on the local acceleration of particles. The time integration was performed with high enough accuracy such that the total energy was conserved to better than 0.1% in all cases. We have also explicitly checked that our results are not compromised by choices of force softening, time-stepping or opening angle criteria in the treecode.

We perform a first experiment in which we merge two identical spherical haloes with zero angular momentum corresponding to a head-on collision with an impact parameter $`b=0`$. The initial relative separation of the two haloes is $`r_{\mathrm{sep}}=210h^1\mathrm{kpc}`$ and are they infalling at a relative speed of $`v_{\mathrm{rel}}=300\mathrm{km}\mathrm{s}^1`$. Relative velocities significantly higher than this value lead to unbound orbits, so we chose this value so as to maximize the amount of flattening from a single merger. Note that the initial separation corresponds to haloes which are overlapping. This does not influence the final result as nearly identical remnants were produced with non-overlapping initial models. The two haloes merge within a time-scale of order the crossing time $`t_{\mathrm{cross}}=\sqrt{r_{\mathrm{vir}}^3/GM_{\mathrm{vir}}}`$ which is approximately equal to $`2`$ Gyr for our initial models. In all simulations the merger remnants were allowed to settle into equilibrium for several crossing times after the merger was complete which was established by monitoring the fluctuations of the density profile of the final system. In addition, the virial ratio $`2T/W`$ varied between 0.1 and 0.3% from $`1.0`$ for the same time-scales confirming that our triaxial remnants are indeed in equilibrium. In all cases the centre of the remnant is identified using the most bound particle, which agrees very well with the centre of mass recursively calculated using smaller spherical regions.

The radial merger produces a prolate axisymmetric system with axial ratios $`a`$:$`b`$:$`c`$ = $`2`$:$`1`$:$`1`$ where $`a`$, $`b`$ and $`c`$ are the long, intermediate and short axis, respectively. For each remnant, we calculate principle axis ratios $`s=b/a`$ and $`q=c/a`$ ($`a>b>c`$), from the eigenvalues of a modified dimensionless inertia tensor (e.g., Dubinski & Carlberg, 1991): $`I_{ij}=x_ix_j/q^2`$, where $`q^2=x^2+(y/s)^2+(z/q)^2`$ is the ellipsoidal coordinate. We use an iterative algorithm starting with a spherical configuration ($`a=b=c`$), and use the results of the previous iteration to define the principle axes of the next iteration. The procedure is iterated until the values of both axial ratios $`b/a`$ and $`c/a`$ have a percentage change of less than $`10^3`$. The dimensionless spin parameter of this remnant, which is defined by

$$\lambda =\frac{J|E|^{\frac{1}{2}}}{GM^{\frac{5}{2}}}$$
(4)
where $`J`$, $`E`$ and $`M`$ are the total angular momentum, binding energy and mass respectively, corresponds to $`\lambda =0`$ as expected in a radial merger.

In a second experiment we start with the two spherical haloes separated again by $`r_{\mathrm{sep}}=210h^1\mathrm{kpc}`$, but we give one halo a transverse velocity equal to the circular velocity of the combined model at that distance, $`v_{\mathrm{rel}}=240\mathrm{km}\mathrm{s}^1`$. The circular orbit merger produces an oblate halo with axial ratios $`a`$:$`b`$:$`c`$ = $`2`$:$`2`$:$`1`$. As we demonstrate later (see Figure 4), haloes with any amount of triaxiality between these values could be created by using unequal-mass mergers or mergers with less or more angular momentum. The spin parameter of the final system is $`\lambda =0.1`$. As expected, with these experiments we have been able to reproduce most of the range in $`\lambda `$ seen in cosmological $`N`$-body simulations (e.g., Warren et al., 1992; Lemson & Kauffmann, 1999). Although the formation histories of halos with a hierarchical universe is significantly more complex than these binary mergers, we might expect that oblate halos have a higher spin parameter and more angular momentum than prolate halos.

Figure 1 shows the density profiles of the initial spherical model and that of the prolate and oblate haloes. The panel on the left shows results using spherically-averaged bins. The density profiles on the right panel were calculated as a function of the ellipsoidal coordinate $`q`$, where $`q^2=x^2+(y/b)^2+(z/c)^2`$. In agreement with other studies (e.g., White, 1978; Villumsen, 1983; Barnes, 1999; Fulton & Barnes, 2001; Boylan-Kolchin & Ma, 2003), remnant central density slopes remain unchanged. In Figure 2 we show the circular velocity curves of the initial model and the final prolate and oblate haloes using again spherically averaged (left panel) and ellipsoidally averaged (right panel) bins. These have all been scaled by the same factor in radius and velocity such that all the curves have the same peak circular velocity. The prolate halo has a nearly identical concentration to the initial spherical model, but the oblate halo is nearly a factor of two less concentrated.

In order to create haloes with higher amounts of flattening it is necessary to merge together the products of the first set of simulations. Merging two prolate haloes along the long axis again at a velocity such that the system is just bound produces a $`a`$:$`b`$:$`c`$ = $`3`$:$`1`$:$`1`$ remnant, while merging two oblate haloes in the plane of the long axis, but counter-rotating, results in a further flattened system with axial ratios $`a`$:$`b`$:$`c`$ = $`3`$:$`3`$:$`1`$ and no net angular momentum, but with particles streaming in opposite directions. Figure 3 shows two contoured images of the oblate halo viewed projected along the short and long axis, respectively.

Finally, we performed a minor merger between two spherical halos with mass ratio equal to $`4:1`$. The most massive halo was identical to the ones used in the previous experiments while the concentration of the smaller one was scaled according to Bullock et al. (2001) for the adopted $`\mathrm{\Lambda }`$CDM cosmology. The orbital angular momentum of the merger was typical of cosmological encounters (Khochfar & Burkert, 2003). The number of particles and gravitational softening used for the less massive halo was such that both halos were resolved with the same mass and force resolution. This is useful in the intepretation of the results regarding the structure of the remnant. Kazantzidis et al. (2004b) examined the evolution of disk galaxy merger remants and found that their shapes evolve in the outer regions for several crossing times. In Figure 4, we demonstrate the stability of this remnant by running the simulation for a further $`10`$ half-mass dynamical times and calculated profiles of axial ratios and triaxiality parameter \[$`T(a^2b^2)/(a^2c^2)`$\] at $`20`$ intermediate steps. The error bars in Figure 4 show the rms variation in the profiles at each radius and demonstrate that the remnant is in equilibrium.

## 3 TIDAL EVOLUTION OF TRIAXIAL SUBHALOES

Numerical investigations of the tidal stripping and mass loss of substructure haloes in CDM models have previously been restricted to spherical systems with isotropic velocity dispersion tensors (e.g., Taffoni et al., 2003; Hayashi et al., 2003). The results of these simulations are also incorporated into semi-analytic models which attempt to follow the tidal evolution of orbiting subhaloes. However, the shapes and velocity dispersion tensors of haloes in cosmological simulations have been studied by many authors (e.g., Warren et al., 1992; Cole & Lacey, 1996; Thomas et al., 1998; Colín, Klypin & Kravtsov, 2000; Bullock, 2002) and they exhibit a significant departure from spherical symmetry and isotropy. It is interesting to investigate if the response of a spherical subhalo to an imposed tidal field is the same as that of a triaxial subhalo.

KMM compared the resilience to tidal forces of two identical self-consistent NFW satellites on the same external tidal field and orbit, but constructed with different velocity dispersion tensors. Identifying the bound mass as function of time, these authors demonstrated that the satellite with a radially anisotropic velocity dispersion (Osipkov, 1979; Merritt, 1985) experiences much more efficient mass loss than its isotropic counterpart. This is due to the fact that particles on more radial orbits spend on average more time at larger radii and are therefore more easily stripped by the external field.

Here we expand on the previous result and show that similar differences are obtained when one considers the tidal evolution of spherical and triaxial haloes as they orbit within a deeper potential. We compare results of the $`3`$:$`1`$ prolate and oblate haloes and the initial isotropic spherical halo used to create these models. The density profile of the spherical model is scaled by a factor of four in mass to preserve the same spherically averaged density profile as the prolate halo. Each halo is then placed on an identical orbit within a static cored isothermal potential such that the ratio of host to satellite halo circular velocities is $`v_{\mathrm{host}}/v_{\mathrm{sat}}`$ = $`10`$:$`1`$ and the apocentric to pericentric radii is $`r_{\mathrm{apo}}/r_{\mathrm{per}}`$ = $`4`$:$`1`$. The alignment of the axes of the triaxial satellite were varied with respect to the external deeper potential and the satellites always begin at the apocentre which is set equal to the virial radius. This allows the satellite-host system to be rescaled to any galactic or cluster system. The imposed potential is reasonably soft since we use a core radius equal to half the pericentric radius. Significantly higher rates of mass loss may be expected in a more cuspy potential, however here we aim to explore the slow tidal loss of mass rather than mass loss due to strong tidal shocks. In a forthcoming paper we will give an analytical model for the rate of mass loss as a function of host potential, orbit and subhalo structure and shape (Kazantzidis et al. , in preparation).

In Figure 5 we perform a comparison of the tidal evolution of each of these three models by identifying the bound mass as a function of time using the publicly available group finder SKID111http://www-hpcc.astro.washington.edu/tools/skid.html (Stadel, 2001). The prolate halo experiences mass loss at a much higher rate than both the spherical and oblate ones. We also found that this result is not significantly sensitive to the alignment between the major or minor axis of the satellite and the orbital plane with the maximum differences in the mass loss being of the order of 25%. In Figures 6 and 7 we show the change in the shape and in the anisotropy parameter $`\beta `$ of our models. We suspect that the prolate halo becomes more spherical due to the tidal removal of particles on radial orbits that are supporting the shape of the system. Future studies will analyse the orbital evolution of particles within the satellite to study these effects in detail. In addition, after $`5`$ Gyr of orbital evolution the velocity dispersion tensor of the prolate halo has significantly changed and now is nearly isotropic.

It is interesting to study the evolution of shape and anisotropy for haloes in a CDM simulation that enter into a deeper potential. This is important when we want to model the orbits of stars in satellite galaxies or the effect of subhaloes on lense models. Here we analyse subhaloes from a high resolution $`\mathrm{\Lambda }`$CDM cluster simulation (Diemand et al., 2004). The initial conditions are generated with the GRAFIC2 package (Bertschinger, 2001). We start with a $`300^3`$ particle cubic grid with a comoving cube size of $`300\mathrm{Mpc}`$ (particle mass $`m_\mathrm{p}=2.6\times 10^{10}h^1\mathrm{M}_{}`$). We trace back and refine a cluster region by a factor of $`12`$ in length and $`1728`$ in mass, so that the mass resolution is $`m_\mathrm{p}=1.5\times 10^7h^1\mathrm{M}_{}`$. The softening length is comoving from the start of the simulation ($`z40`$) to $`z=9`$. From $`z=9`$ until present, we use a physical softening length of $`1.05\times 10^3r_{\mathrm{vir},\mathrm{parent}}`$. At $`z=0`$ the refined cluster contains $`1.4\times 10^7`$ particles (R12) and we select a relatively massive subhalo which is resolved with more than $`2\times 10^5`$ particles.

Figure 8 shows the axial ratios $`b/a`$ and $`c/a`$ (left panel) and the anisotropy parameter $`\beta `$ (right panel) as a function of radius for this object. The thin lines correspond to the subhalo just before entering its host, whereas the thick lines correspond to the third pericentre passage of the subhalo. Note that the radius of the subhalo is given in units of $`r_\mathrm{t}`$, where $`r_\mathrm{t}`$ is its tidal radius after the third pericentre passage. Figure 8 confirms that subhaloes become more spherical and their velocity distribution more isotropic owing to tidal effects. We then compare the mean axial ratios of field haloes versus subhaloes in the same simulation. On average we find that subhaloes have axial ratios that are 30% larger (more spherical) than field haloes.

## 4 CONCLUSIONS

We have presented a method for constructing cuspy axisymmetric and triaxial $`N`$-body haloes based on merging isotropic equilibrium spherical haloes with varying amounts of angular momentum. In particular, we found that radial mergers produce prolate systems, while mergers on circular orbits produce oblate systems. This technique has the benefit that it is based on the way in which haloes obtain their triaxiality in cosmological simulations, therefore the anisotropy distributions and angular momentum of haloes are well motivated.

Mergers between similar equilibrium spherical haloes at high resolution show that the resulting halo has the same density profile, independent of the angular momentum of the merger. This has two implications: (i) the density profile of the triaxial halo can be set by the choice of the density profiles of the progenitor haloes and (ii) similar mass mergers can not dramatically re-arrange the central density structure of DM haloes, albeit oblate haloes (high angular momentum merger remnants) have concentrations up to a factor of two lower than prolate haloes (low angular momentum merger remnants). This may explain the entire scatter in the distribution of halo concentrations from cosmological simulations (Bullock et al., 2001; Eke et al., 2001; Wechsler et al., 2002). It may also be the case that the least concentrated haloes host the low surface brightness discs since these galaxies may form in high angular momentum oblate haloes (e.g., O’Nell et al., 1997). These galaxies generally require low values of the concentration when fit to cuspy halo models (e.g., McGaugh & de Blok, 1998; van den Bosch et al., 2000; de Blok et al., 2001; Swaters et al., 2003).

As an application we considered the tidal evolution of triaxial subhaloes orbiting within deeper potentials. Haloes with identical spherically averaged density profiles and on identical orbits evolve self-similarly with time. Spherical haloes with isotropic velocity dispersion tensors suffer significantly less mass loss than radially anisotropic prolate ones. This is likely due to the fact that particles on radial orbits are easily stripped which also results in subhaloes becoming more spherical. In general, the subhaloes become more spherical and their velocity distribution more isotropic owing to tidal effects. This result has been confirmed by investigating the response to tides of both isolated satellite haloes orbiting within a static host potential and substructure haloes in the time dependent cosmological tidal field. We find that subhaloes in a cosmological simulation of a cluster are on average 30% rounder than their field counterparts. Galaxy subhaloes should be even closer to spherical since they spend longer being reshaped by the host potential.

Our more realistic modelling of DM haloes is important for studies of the weak and strong lensing statistics attempting to distinguish between competing cosmological models, the structural evolution and mass loss from substructure, the formation of tidal streams or the sizes of satellite haloes in galaxies and clusters. These results may also be important to incorporate within semi-analytic models that attempt to model the distribution of satellites in DM haloes. CDM haloes are generally radially anisotropic therefore they lose mass and are disrupted more quickly than the isotropic systems that have been generically considered in previous studies. This may explain why semi-analytic models predict far more substructures than found in high resolution numerical studies (e.g., Zentner & Bullock, 2003, James Taylor private communication).

## Acknowledgments

It is a pleasure to thank the anonymous referee for constructive comments on the manuscript and Victor Debattista and Lucio Mayer for useful discussions. BM thanks Priyamvada Natarajan for organising the Yale Cosmology Workshop “The Shapes of Galaxies and Their Dark Matter Haloes” (2001) which motivated some of this work and also apologises for not writing up the conference proceedings. The numerical simulations were carried out on the zBox (http://www-theorie.physik.unizh.ch/ stadel/). Initial conditions for the cosmological simulations were generated at the Swiss Center for Scientific Computing (SCSC) at Manno.


text:
# On the reconstruction of diffusions from first-exit time distributions

## 1 Introduction

Consider a diffusive stochastic process $`X_t`$, with an initial condition $`X_0`$. Once $`X_t`$ reaches a certain threshold value $`X^{}`$, the process is stopped and restarted at $`X_0`$. The distribution of times for $`X_t`$ to first reach $`X^{}`$ across the ensemble of measurements, the first-exit time distribution, arises in many physical applications . Mean first-exit times are readily calculated from explicit forms for the dynamical rules controlling the process $`X_t`$. However, in many applications, the underlying physics is unknown, complex, and/or very difficult to model. In such instances, the question arises as to whether one can reconstruct the form of the governing equations or equation parameters that best fit measured data (such as first-exit time distributions). In this paper, we analyse a general scalar diffusive process to formally determine how much of the drift and diffusion functions can be determined from first-exit time data.

Many first passage time problems can be recast as inverse problems. We describe two physical realizations that motivate this work: transmembrane voltage spike frequencies in neurons , and first passage times over a molecular barrier in the rupture of chemical bonds . The transmembrane voltage of a nerve cell is routinely measured and typically exhibits repeated, sharp spikes. A class of reduced models for the transmembrane voltage $`V`$ is described by

$$\frac{dV_t}{dt}=I(V_t)+g(V_t)\eta _t,$$
(1)
where $`I(V_t)`$ and $`g(V_t)`$ are typically polynomial functions of $`V`$, and $`\eta _t`$ is a noise arising from interactions with the many other connected neurons. This Langevin-type equation is to be appropriately interpreted with Stratonovich calculus . An additional, instantaneous nonlinearity that gives rise to a spike and voltage reset is implicitly included by imposing a threshold voltage $`V^{}`$. This nonlinearity arises ultimately from the nonlinear dynamics of ion channels that span the cell membrane . Upon reaching $`V^{}`$, the system instantaneously spikes and resets to the value $`V_0`$. Various forms of $`I(V_t)`$ and $`g(V_t)`$ have found wide use. For example, when the noise $`\eta _t`$ arising from other connected nerve cells changes the nerve cell’s membrane conductivity, $`g(V_t)V_t+\text{const.}`$ Both linear (linear integrate-and-fire) and quadratic (quadratic integrate-and-fire) forms are used to model $`I(V_t)`$ , although it is often unclear which model of $`I(V_t)`$ is most appropriate . Instead of solving for the spike times for presupposed $`I(V_t)`$ and $`g(V_t)`$ , our goal is to infer as much as possible about the nonlinear current-voltage relationship $`I(V_t)`$ and coupling to noisy inputs $`g(V_t)`$, given a measured spike time distribution .

Another application where first exit times can be used to reconstruct models is in the single molecule measurements of macromolecular detachment, or bond rupturing. An example of this rapidly developing field is the accelerated rupturing of bonds upon application of a dynamic load between two macromolecules, such as the ligand-receptor complex biotin-streptavidin . The mean times to bond breaking have been measured  and modelled . This “Kramers problem” , can also be studied as an inverse problem. Bond rupture can be modelled as the irreversible passage of the relevant molecular coordinate $`X`$ over a barrier from a metastable state. Thermal bombardments agitate both the molecular distance coordinate directly and the interaction potential. In the overdamped, or Markovian limit, the distance coordinate obeys the Langevin equation

$$\frac{dX_t}{dt}=U^{}(X_t)+g(X_t)\eta _t,$$
(2)
where $`U(X_t)`$ is the effective interaction potential of the rupturing bond, and $`g(X_t)`$ represents a force that arises from a fluctuating potential driven by the white diffusion noise $`\eta _t`$ (again, the appropriate convention is the Stratonovich difference rule). The forward problem of bond rupture, or detachment, using simple fixed potentials $`U(X_t)`$ , and with explicit forms for fluctuating potentials , have been well-studied.

Both problems discussed above (see Fig. 1), as well as many others, can be modeled using a scalar stochastic process where the first-exit time distribution obeys the backward Kolmogorov equation. We show in the following sections, that the measured first-exit time distribution (given a fixed initial condition) along with knowledge of the drift and diffusion parameters in part of the domain, is sufficient to uniquely determine one of them in the entire domain provided the other one is known. Our analysis gives a formal procedure for reconstructing drift or diffusion of the stochastic process. It is based on transforming the equation into a Schrödinger type equation and on analytical (more precisely meromorphic) continuation of measured data to obtain the spectral data of the Schrödinger equation. Although uniqueness is demonstrated in the aforementioned framework, the necessary analytical continuation in the reconstruction shows that the inverse problem is severely ill-posed .

## 2 Stochastic Equations and Main Result

Consider the general stochastic process

$$\begin{array}{c}dX_t=\mu (X_t)dt+\sigma (X_t)dW_t,\hfill \end{array}$$
(3)
where $`dW_t=\eta _tdt`$ is the Wiener measure, and the Itô interpretation is used . The drift term, $`\mu (X)`$, is assumed to have sub-linear growth at infinity, and is of class $`C^1()`$. We also assume that the variance, $`\sigma (X)`$, is uniformly bounded from above and below by positive constants, and is of class $`C^2()`$. These assumptions are sufficient to ensure that the process (3) admits a unique solution .

In the application to spiking action potentials in neurons, $`X_tV_t`$ is the transmembrane voltage at time $`t`$. The potential immediately after reset is denoted $`X_0V_0`$. The corresponding drift and diffusion terms are defined as $`\mu (X)I(X)+\frac{1}{2}g(X)g^{}(X)`$ and $`\sigma (X)g(X)`$, respectively, since (3) is in the Itô form and (1) should be interpreted in the Stratonovich sense . For a more detailed derivation of the the diffusion process in the context of action potential modeling, we refer the reader to number of standard works . In the application to chemical bond rupturing, $`X`$ denotes the reaction coordinate along an effective molecular energy landscape (cf. Fig. 1(b)), $`X_0`$ is an initial bond displacement, and we identify the drift in (3) with $`U^{}(X)+\frac{1}{2}g(X)g^{}(X)`$.

Let us now consider the interval $`\mathrm{\Delta }=(X_{},X^{})`$ and the first exit time $`t^{}`$ of the process $`X_t`$; i.e., $`t^{}`$ is the first time of exit of $`X_t`$ out of $`(X_{},X^{})`$ . First exit time measurements can then be used to reconstruct properties of the diffusion process. For instance, it is known that knowledge of $`P_X(X_t^{}=X^{})`$ (the probability that $`X_t`$ exits $`\mathrm{\Delta }`$ at $`X^{}`$ conditional on it starting at $`X\mathrm{\Delta }`$) and $`𝐄_X[t^{}]`$ (the expected time it takes for $`X_t`$ to exit $`\mathrm{\Delta }`$ starting at $`X`$), uniquely determines $`\mu `$ and $`\sigma `$. We refer to Breiman  for details on the so-called natural time scale and speed measure of diffusion processes.

In the applications presented here, the situation is different: We have access to the full distribution of the first-exit times, but only for processes starting at $`X=X_0`$ (the reset potential or the initial relative molecular coordinate). The distribution of the exit times can be conveniently written as the solution to a partial differential equation (the backward Kolmogorov equation). Let us denote by $`w(X;\lambda )`$ the Laplace transform (from the time variable $`t`$ to the dual variable $`\lambda `$) of the probability density $`\stackrel{~}{w}(X;t)`$ of the first-exit times $`t`$ from $`X\mathrm{\Delta }`$ of the diffusion process starting at $`X=X_0`$. For simplicity, we will assume reflecting boundary conditions for the diffusion at $`X=X_{}`$, and consider termination of the process only at $`X^{}`$. Following , we find that $`w(X;\lambda )`$ solves the following backward Kolmogorov equation,

$$\begin{array}{c}\frac{\sigma ^2(X)}{2}w^{\prime \prime }(X;\lambda )+\mu (X)w^{}(X;\lambda )=\lambda w(X;\lambda ),X_{}<X<X^{},\lambda >0,\hfill \\ w(X^{};\lambda )=1,\frac{w(X_{};\lambda )}{X}=0,\hfill \end{array}$$
(4)
where $`w(X_{};\lambda )/X=0`$ corresponds to vanishing probability flux at $`X=X_{}`$. Although we restrict our analysis to a reflecting lower bound, depending on the physical application, other boundary conditions can be imposed at an arbitrary lower bound. The boundary condition $`w(X^{};\lambda )=1`$ reflects the fact that $`\stackrel{~}{w}(X_0=X^{},t)=\delta (t)`$. Since $`\stackrel{~}{w}(X_0;t)`$ is directly measured, $`w(X_0;\lambda )`$ is exactly determined for all $`\lambda >0`$. Using an affine change of variables, $`(XX_{})/(X^{}X_{})x`$, we rescale the domain to $`x[0,1]`$ and recast Eq. (4) as

$$\begin{array}{c}\frac{\sigma ^2(x)}{2}w^{\prime \prime }(x;\lambda )+\mu (x)w^{}(x;\lambda )=\lambda w(x;\lambda ),0<x<1,\lambda >0,\hfill \\ w(1;\lambda )=1,\frac{w(0;\lambda )}{x}=0,\hfill \end{array}$$
(5)
where $`x`$ is now the initial position of the process. We still use the notation $`\sigma `$ and $`\mu `$ for the diffusion and drift terms in the rescaled variables. For example, in the neuron spiking problem, the voltage threshold is at $`1`$, the inhibitory threshold (the lower bound for the voltage) is at $`0`$, and the reset voltage that restarts the process is at $`x=x_0`$. This problem admits a unique solution provided $`\sigma `$ remains uniformly bounded from above and below and $`\mu `$ remains bounded.

In many instances $`\mu (x)`$ and $`\sigma (x)`$ are known from a well-accepted model, and the function $`w(x;\lambda )`$ can be computed directly. Now consider the case where $`\mu (x)`$ and $`\sigma (x)`$ are unknown, but $`w(x=x_0;\lambda )`$ is measured at a single initial position $`x_0=(X_0X_{})/(X^{}X_{})`$. What can one reconstruct (with respect to $`\mu (x)`$ and $`\sigma (x)`$) from a knowledge of $`w(x_0;\lambda )`$ for all $`\lambda >0`$? In answering this question, we arrive at the following:

Let $`0<x_b<1`$ be the point defined in (18) below and let $`\mathrm{\Lambda }(0,1)`$ be a countable set of Lebesgue measure $`0`$ in $`(0,1)`$, defined in (23) below. We obtained the following related results:

###### Theorem 1

Assume that $`\mu (x)`$ and $`\sigma (x)`$ are known on $`(0,x_b)`$ and that the measurement point $`x_0(0,1)\backslash \mathrm{\Lambda }`$. If, in addition, either $`\mu (x)`$, $`\sigma (x)`$, or $`\mu (x)/\sigma (x)`$ is known entirely on $`[0,1]`$, then both $`\mu (x)`$ and $`\sigma (x)`$ are uniquely determined on $`[0,1]`$ by measurement of $`w(x_0;\lambda )`$, for all $`\lambda >0`$.

###### Theorem 2

Assume that $`\mu (x)`$ and $`\sigma (x)`$ are known on $`(0,x_0)`$ and that the measurement point $`x_0x_b`$. If, in addition, either $`\mu (x)`$, $`\sigma (x)`$, or $`\mu (x)/\sigma (x)`$ is known on $`[0,1]`$, then both $`\mu (x)`$ and $`\sigma (x)`$ are uniquely determined on $`[0,1]`$ by measurement of $`w(x_0;\lambda )`$, for all $`\lambda >0`$.

Both theorems say that if either $`\mu `$ or $`\sigma `$ is known on $`[0,1]`$, while the other is known only on half of the domain, we can reconstruct this other function also on the entire domain. Since $`x_b`$ is implicitly defined by $`\sigma `$ (cf. Eq. 18), it is an a priori constraint on the size of the domain that one can reconstruct from the measurements $`w(x_0;\lambda )`$. In Theorem 1, the measurements are made at an almost arbitrary point $`x_0(0,1)`$. We have to remove a countable set of points, which is also defined implicitly as it depends on $`\sigma `$ and $`\mu `$. This constraint is no longer necessary in Theorem 2, when $`x_0(x_b,1)`$.

The conclusions of both results hold true if the domains $`(0,x_b)`$ and $`(0,x_0)`$ of a priori knowledge of the coefficients are replaced by $`(x_b,1)`$ and $`(x_0,1)`$, respectively with $`x_0x_b`$ in the latter case. We do not explicitly consider these cases. The main steps of the proof are the following. First we use a transformation similar to the Liouville transformation (used to obtain the canonical form of Sturm-Liouville problems ) to map Eq. (5) into a Schrödinger equation (Section 3). Next, we show that a measurement of $`w(x_0;\lambda )`$, for all $`\lambda >0`$, is sufficient to obtain a single eigenvalue spectrum of the Sturm-Liouville problem (when $`x_0\mathrm{\Lambda }`$ in the setting of Theorem 1). This allows us to use results from Gelfand and Levitan  to verify uniqueness of the reconstruction (Section 4) and explicitly extract $`\mu (x)`$ or $`\sigma (x)`$ (Section 5). The above theorems require to know the coefficients on half of the domain before one of them can be reconstructed on the other half. An explicit asymptotic expansion presented in section 5 shows that this a priori knowledge is indeed necessary to obtain a unique reconstruction. In section 6 we present a result showing that the reconstruction of one of the coefficients $`\mu (x)`$ or $`\sigma (x)`$ is unique over the whole domain $`(0,1)`$ provided that they are known on an arbitrarily small domain $`(0,\eta )`$ and provided that a sequence of measurements $`w(x_m;\lambda )`$ is available for $`\lambda >0`$ and $`x_m`$, $`m`$, judiciously chosen. In this result, a minimum of $`M`$ measurements is required to reconstruct $`\mu `$ or $`\sigma `$ provided that they are known a priori on an interval close to $`x=0`$ of size $`2^M`$. Finally section 7 considers possible implementations of the method.

## 3 Mapping to Schrödinger form

To recast Eq. (5) as a Schrödinger equation, we effect the following change of variables

$$w(x)=(fu)h(x).$$
(6)
Thus,

$$w^{}(x)h^{}(fu)^{}h,w^{\prime \prime }(x)h^{\prime \prime }(fu)^{}h+(h^{})^2(fu)^{\prime \prime }h,$$
(7)
and Eq. (5) becomes

$$\left[\frac{\sigma ^2}{2}(h^{})^2\right](fu)^{\prime \prime }h+\left[\frac{\sigma ^2}{2}h^{\prime \prime }+\mu h^{}\right](fu)^{}h=\lambda (fu)h.$$
(8)
To normalise the second-order term, we impose the constraint

$$\begin{array}{c}h^{}(x)=\frac{\sqrt{2}}{\sigma (x)},0<x<1,\hfill \\ h(0)=0.\hfill \end{array}$$
(9)
Since $`\sigma `$ is bounded from above and below by a positive constant, $`xh(x)`$ is a diffeomorphism. Upon introducing the change of variables

$$y=h(x),$$
(10)
the points $`x=0`$, $`x=x_0`$, and $`x=x_b`$ before the change of variables become $`y=h(0)=0`$, $`y_0h(x_0)`$, and $`y_bh(x_b)`$, respectively. If we define a new drift

$$\nu (y)=\left(\frac{\sigma ^2}{2}h^{\prime \prime }+\mu h^{}\right)(h^1(y)),0yy^{}h(1),$$
(11)
Eq. (8) can be succinctly expressed as

$$(fu)^{\prime \prime }(y)+\nu (y)(fu)^{}(y)=\lambda (fu)(y).$$
(12)
To remove the drift term, we force $`f`$ to be the solution to

$$\begin{array}{c}2f^{}(y)+\nu (y)f(y)=0,0<y<y^{},\hfill \\ f(y_0)=1.\hfill \end{array}$$
(13)
Since $`f(y)>0`$ for $`y=y_0`$, $`f(y)>0`$ for all $`y[0,y^{}]`$. We finally obtain

$$\begin{array}{c}u^{\prime \prime }(y)q(y)u(y)=\lambda u(y),0<y<y^{}\hfill \\ u(y^{})=\frac{1}{f(y^{})}>0,\frac{u(0;\lambda )}{y}=0,\hfill \end{array}$$
(14)
where the potential is defined by

$$q(y)=\frac{\nu f^{}+f^{\prime \prime }}{f}(y).$$
(15)
By assumption, $`\sigma (x)`$ and $`\mu (x)`$ are known for $`0xb`$. This implies that $`\nu (y)`$ and $`f(y)`$, and hence $`q(y)`$, are known for $`0yy_bh(x_b)`$, although $`y^{}h(1)`$ itself is generally not known (unless $`\sigma (x)`$ is also known on $`[0,1]`$). This undetermined domain size is found in the following section.

## 4 Uniqueness of the potential reconstruction

We now consider the reconstruction of the potential $`q(y)`$ in the equation

$$\begin{array}{c}u^{\prime \prime }(y;\lambda )q(y)u(y;\lambda )=\lambda u(y;\lambda ),0<y<y^{}h(1)\hfill \\ u(y^{})=\frac{1}{f(y^{})}>0,\frac{u(0;\lambda )}{y}=0.\hfill \end{array}$$
(16)
The assumption is that $`q(y)`$ is known on $`[0,y_b]`$ as described earlier. However $`y^{}=h(1)`$ depends on $`\sigma (x)`$ throughout the domain and is unknown, as is $`q(y)`$ on $`(0,y^{})`$. From measurements, we have $`u(y_0;\lambda )`$ for all $`\lambda >0`$ since $`w(x_0;\lambda )=u(y_0;\lambda )`$ and $`f(y_0)=1`$.

Uniqueness of the reconstruction of $`q(y)`$ on $`(y_b,y^{})`$ is based on the fact that $`\lambda u(y;\lambda )`$ is analytic with poles corresponding to values of $`\lambda `$ that are eigenvalues for specific spectral equations. It is known that two spectral equations are necessary to reconstruct the potential $`q(y)`$ from spectral data . However, we will see that our problem allows extraction of at most only one spectrum. Nonetheless, if the “potential” $`q(y)`$ is known for half of the entire domain, one eigenvalue spectrum suffices for the complete and unique reconstruction of $`q(y)`$ .

Since $`q(y)`$ is known for $`0yh(x_b)`$, we require

$$h(x_b)=\frac{y^{}}{2}\frac{h(1)}{2}$$
(17)
in order for $`q(y)`$ to be completely reconstructed on $`[0,y^{}]`$ Therefore, the original functions $`\mu (x)`$ and $`\sigma (x)`$ must at least be known for $`0xx_b`$, where

$$x_b=h^1\left(\frac{y^{}}{2}\right)=h^1\left(\frac{h(1)}{2}\right).$$
(18)
The full problem on $`[0,y^{}]`$ becomes

$$\begin{array}{c}u^{\prime \prime }(y)q(y)u(y)=\lambda u(y),0<y<y^{}h(1)\hfill \\ u(y^{})=\frac{1}{f(y^{})}0,\frac{u(0;\lambda )}{y}=0.\hfill \\ u(y_0;\lambda ),\text{ known.}\hfill \end{array}$$
(19)
Now, consider the associated eigenvalue problem

$$\begin{array}{c}\varphi _k^{\prime \prime }(y)q(y)\varphi _k(y)=\lambda _k\varphi _k(y),0<y<y^{},k1\hfill \\ \varphi _k^{}(0)=\varphi _k(y^{})=0.\hfill \end{array}$$
(20)
If the $`\lambda _k`$ in (20) can be determined from our measured data $`u(y_0;\lambda )`$, then $`q(y)`$ on $`[0,y^{}]`$ can be determined upon using classic theories described in . To find the eigenvalues to (20), as well as the transformed threshold $`y^{}`$, we decompose the Green’s function for (19) in terms of the eigenfunctions of (20), and obtain

$$u(y;\lambda )=\underset{k}{}\left[\frac{\varphi _k(y)\varphi _k^{}(y^{})u(y^{})}{\lambda +\lambda _k}\right].$$
(21)
First, consider the setting of Theorem 1. Since $`u(y^{})=1/f(y^{})>0`$ and $`|\varphi _k^{}(y^{})|>0`$ (since $`\varphi _k`$ is a non-trivial solution to (20)), $`|u(y_0;\lambda )|\mathrm{}`$ as $`\lambda \lambda _k`$, provided $`\varphi _k(y_0)0`$. The latter constraint is not satisfied for all choices of $`y_0`$. Let us define

$$\mathrm{\Lambda }_y=\{y(0,y^{}),\text{ there exists }k1\text{ such that }\varphi _k(y)=0\}.$$
(22)
Since the number of zeros of solutions to the Sturm-Liouville problem (20) is finite, the set $`\mathrm{\Lambda }_y`$ is countable and therefore of Lebesgue measure $`0`$ on $`(0,y^{})`$. As long as our initial measurement position $`y_0`$ is not part of the above countable subset of values where the eigenfunctions of (20) have a node, we can in principle obtain one spectrum by analytically continuing the data $`u(y_0;\lambda )`$ in $`\lambda `$ and finding the poles. The constraint on $`y_0`$ also amounts to saying that the eigenfunctions $`\varphi _k`$ cannot be restricted to being eigenvectors of the same problem posed on the domain $`(0,y_0)`$ instead of $`(0,y^{})`$. In the primitive variables the countable set of forbidden measurement points is thus defined as

$$\mathrm{\Lambda }=h^1(\mathrm{\Lambda }_y)=\{x(0,1),h(x)\mathrm{\Lambda }_y\}.$$
(23)
Now consider Theorem 2. We can extend $`u(y;\lambda )`$ by evenness on $`(y^{},y^{})`$ and then by periodicity on $``$. The resulting function is then continuous in $`y`$ and we deduce from (21) that

$$\frac{u}{y}(y;\lambda )=\underset{k}{}\left[\frac{\varphi _k^{}(y)\varphi _k^{}(y^{})u(y^{})}{\lambda +\lambda _k}\right].$$
(24)
Now since $`\mu (x)`$ and $`\sigma (x)`$ are known on $`(0,x_0)`$ then so is $`q(y)`$ on $`(0,y_0)`$. Since $`u(y_0,\lambda )`$ is known, we can solve (14) for $`u(y)`$ on $`(0,y_0)`$ and thus obtain the value of $`\frac{u}{y}(y_0;\lambda )`$. Since either $`|\varphi _k(y_0)|>0`$ or $`|\varphi _k^{}(y_0)|>0`$ (for otherwise $`\varphi _k(y)0`$), we deduce that $`|u(y_0;\lambda )|+|\frac{u}{y}(y_0;\lambda )|\mathrm{}`$ as $`\lambda \lambda _k`$. Under the hypotheses of Theorem 1 and Theorem 2, we thus find that one can reconstruct the spectrum $`\{\lambda _k\}_{k1}`$ in (20) from the measurements.

Finally, we must determine $`y^{}`$. Asymptotic expansions of the eigenvalues for large $`k`$, found using the WKB approximation  have the form:

$$\lambda _k=\left[\left(k\frac{1}{2}\right)\frac{\pi }{y^{}}\right]^2+\frac{1}{y^{}}_0^y^{}q(y)𝑑y+a_k,$$
(25)
where the sequence $`a_kl^2`$, i.e., $`a_k=(_{k=1}^{\mathrm{}}a_k^2)^{1/2}<\mathrm{}`$. Equation (25) implies that the threshold value $`y^{}`$ can be evaluated using

$$y^{}=\underset{k\mathrm{}}{lim}\frac{k\pi }{\sqrt{\lambda _k}}.$$
(26)
Thus, the domain size is reconstructed. The classical theories found in  can then be directly applied to uniquely reconstruct the remaining, unknown part of $`q(y)`$. Moreover, efficient numerical algorithms have been developed for reconstructing $`q(y)`$ , and ultimately $`\mu `$ and $`\sigma `$ from the spectral data. However, the reconstruction of these spectral data from measurements is based on analytic continuation. This implies that the inverse problem is severely ill-posed  and that noise in the data is quite strongly amplified in the reconstruction; see section 7. Note that a similar approach was used in  in the context of the one-dimensional heat equation.

## 5 Reconstruction of drift or diffusion term

We showed in the preceding section that $`q(y)`$ and $`y^{}`$ are uniquely determined from the measurements provided that one set of hypotheses stated in Theorems 1-3 (cf. Section 6) are verified. Combining equations (13) and (15), we find

$$\nu ^{}(y)+\frac{\nu ^2}{2}=2q(y),0<y<y^{}.$$
(27)
Since both $`\mu `$ and $`\sigma `$ are known on $`[0,x_b]`$,

$$\nu (y_b)=\left(\frac{\sqrt{2}\mu (x_b)}{\sigma (x_b)}\frac{\sigma ^{}(x_b)}{\sqrt{2}}\right).$$
(28)
The Riccati equation (27), along with the boundary condition (28) admits a unique solution for $`\nu (y)`$.

Even if $`\nu (0yy^{})`$ is known uniquely, one cannot extract both $`\mu (x)`$ and $`\sigma (x)`$ independently. However, if either $`\mu `$ or $`\sigma `$ is known, then the other one can uniquely be reconstructed from the knowledge of $`\nu (y)`$ on $`(0,y^{})`$. The simplest case occurs when only $`\sigma `$ is known, for then, so is $`h(x)`$, and the reconstruction of $`\mu `$ follows from the definition (11) of $`\nu `$. On the other hand, we can deduce from (9) and (11) that

$$\left(\frac{\sigma ^{}}{2}\frac{\mu }{\sigma }\right)(x)=\frac{1}{\sqrt{2}}\nu \left(_0^x\frac{\sqrt{2}}{\sigma (x^{})}𝑑x^{}\right).$$
(29)
Let us assume that $`\mu /\sigma `$ can be expressed as a known (smooth) function of $`\sigma `$ and $`x`$. This includes the cases when $`\mu `$ is known or when $`\mu /\sigma `$ is known. Since $`\sigma `$ is a priori known to be uniformly bounded from below, equation (29) provides a first-order integro-differential equation that can then be solved for $`\sigma (x)`$ on $`(0,1)`$. The reconstruction is uniquely defined, concluding the proof of the main theorem.

The reconstruction of $`q(y)`$ on $`y(0,y^{})`$ clearly depends on the choice of $`q(y)`$ on $`[0,y_b]`$ since it would change the spectrum through Eq. (20). Therefore, a physically incorrect model of $`q`$ on $`[0,y_b]`$, contaminates the reconstruction of $`\mu (x),\sigma (x)`$ on $`x(0,1)`$, leading to physically incorrect, albeit unique, results. To be more specific, consider the simple case $`\sigma ^22`$, whence $`h(x)=x`$, $`y^{}=1`$, and $`\nu (x)=\mu (x)`$. Provided $`|\mu |^2|\mu ^{}|`$, Eq. 15 implies $`q(x)\mu ^{}(x)/21`$. We then perform an asymptotic expansion of the eigenfunctions $`\varphi _k=\varphi _{k0}+\varphi _{k1}`$ and eigenvalues $`\lambda _k=\lambda _{k0}+\lambda _{k1}`$, where $`\varphi _{k1}1`$ and $`\lambda _{k1}1`$ (see  for instance), and obtain

$$\lambda _{k0}=\left(\pi \left(k\frac{1}{2}\right)\right)^2,\text{ and }\varphi _{k0}(x)=\mathrm{cos}\pi \left(k\frac{1}{2}\right)x,$$

and using the Fredholm alternative in the perturbation

$$\varphi _{k1}^{\prime \prime }+q\varphi _{k0}=\lambda _{k1}\varphi _{k0}+\lambda _{k0}\varphi _{k1},$$

we obtain

$$\lambda _{k1}=\frac{{\displaystyle _0^1}q(x)\varphi _{k0}^2(x)𝑑x}{{\displaystyle _0^1}\varphi _{k0}^2(x)𝑑x}.$$

Using the relation $`\mathrm{cos}2\theta =2\mathrm{cos}^2\theta 1`$, we deduce from the above equation that the knowledge of $`\lambda _{k1}`$ allows us to obtain

$$q_{2m1}=\frac{1}{2}_0^1q(x)\mathrm{cos}(2m1)\pi x𝑑x.$$

Therefore the knowledge of one spectrum $`\lambda _k`$ gives us half of the cosine Fourier transform of $`q(x)`$, hence of $`\mu (x)`$. The other half needs to be known, for instance by imposing that we know $`q(x)`$, or equivalently $`\mu (x)`$, on $`(0,1/2)`$.

Although we have shown what parts of $`\mu (x)`$ and $`\sigma (x)`$ can be uniquely reconstructed from a single first-exit time distribution, in practice, finding the eigenvalues $`\lambda _k`$ from data to the precision required for accurate reconstruction is difficult.

## 6 Reconstruction with multiple measurements

We have seen in the preceding section that either $`\mu `$ or $`\sigma `$ could be reconstructed from the knowledge of $`q`$ provided the other one is known. In this section we obtain sufficient conditions on the measurements so that $`q(x)`$ can uniquely be reconstructed on $`(\epsilon ,1)`$ for $`\epsilon >0`$. We show that it is sufficient to measure $`u(y_m;\lambda )`$ at a countable sequence of points $`y_m`$ with accumulation point at $`y=0`$. These points $`y_m`$ correspond to points $`x_m`$ via the transform $`x_m=h^1(y_m)`$. They are constructed as follows.

We choose $`y_0`$ such that $`y_0y^{}/2`$ and then by induction the points $`y_m`$ ($`m1`$) such that $`y_m(y_{m1}/2,y_{m1})`$. We assume that $`y_m0`$ as $`m\mathrm{}`$. We denote by $`𝒴`$ the set of sequences $`\{y_m\}_m=(y_0,y_1,\mathrm{})`$ satisfying the above constraints. We can then show the following result

###### Theorem 3

Let us assume that $`qL^2(0,1)`$, that there is a small interval $`[0,\epsilon ]`$ where $`q(y)`$ is constant and that there is a known index $`M`$ such that $`y_M<\epsilon `$. The parameters $`\epsilon `$ and $`y_M`$ need not be known. Then the measurements at $`u(y_m;\lambda )`$, $`1mM+1`$, for $`\{y_m\}𝒴`$ and $`\lambda >0`$ uniquely determine $`q(y)`$ on $`(0,y^{})`$ as well as the point $`\{y_m\}_m`$ and $`y^{}`$.

Proof. By hypothesis, we know $`M`$ such that $`y_M<\epsilon `$. Consider the value of $`u(y_M;\lambda )0`$. If $`u(y_M;\lambda )=0`$, the condition at $`y=0`$ would imply that $`u(y;\lambda )0`$ on $`(0,y_M)`$ and, hence on $`(y_M,y^{})`$, since then $`u(y_M;\lambda )=\frac{u}{y}(y_M;\lambda )=0`$. This however contradicts the fact that $`u(y^{})=1`$. Upon dividing the solution $`u(y;\lambda )`$ by $`u(y_M;\lambda )0`$ we deduce that $`v_M(y;\lambda )=u(y;\lambda )/u(y_M;\lambda )`$ solves

$$\begin{array}{c}v_M^{\prime \prime }(y;\lambda )q(y)v_M(y;\lambda )=\lambda v_M(y;\lambda ),0<y<y_M\hfill \\ v_M(y_M;\lambda )=1,\frac{v_M(0;\lambda )}{y}=0.\hfill \end{array}$$
(30)
We then know from the proof of Theorem 2 and the asymptotic expansion (25) that the measurement of $`v(y_{M+1})`$ uniquely determines $`y_M`$ and $`q(y)`$ on $`(0,y_M)`$ since the latter is constant.

The rest of the proof follows by induction. Let us assume that $`q(y)`$ is known on $`(0,y_{m+1})`$ for $`m+1M`$. Then $`u(y_m;\lambda )0`$ and $`v_m(y;\lambda )=u(y;\lambda )/u(y_m;\lambda )`$ satisfies the same equation (30) with $`y_M`$ replaced by $`y_m`$. Upon using Theorem 2 we deduce that the knowledge of $`v_m(y_{m+1};\lambda )`$ for $`\lambda >0`$ uniquely determines $`y_m`$ and $`q(x)`$ on $`(0,y_m)`$. This allows us to reconstruct $`y_0`$ and $`q(y)`$ on $`(0,y_0)`$ by induction. Applying Theorem 2 on $`(0,y^{})`$ one last time concludes the proof of the theorem.

The same type of proof could be used with other local a priori information on the behavior of $`q(y)`$ in the vicinity of $`y=0`$. For instance, we can show that the reconstruction is feasible provided that $`q(x)`$ is known on an arbitrarily small interval $`(0,\epsilon )`$ and provided that there exists an admissible sequence $`(y_0,y_1,\mathrm{},y_M)`$ of measurements with $`y_M\epsilon `$.

In both cases the reconstruction of $`q(y)`$ and of the measurement points $`y_m`$ allows us to obtain some information about the original drift and diffusion terms $`\mu (x)`$ and $`\sigma (x)`$. We can for instance prove the following result.

###### Theorem 4

Assume that $`\mu (x)`$ and $`\sigma (x)`$ are known on $`(0,\eta )`$ for some $`\eta >0`$ and that $`w(x_m;\lambda )`$ is measured for $`\lambda >0`$ at $`x_m`$, $`1mM`$, where the constraints on $`\{x_m\}_m`$ are that $`x_M\eta `$ and $`\{y_m\}_m=\{h(x_m)\}_m𝒴`$. If, in addition, either $`\mu (x)`$, $`\sigma (x)`$, or $`\mu (x)/\sigma (x)`$ is known on $`[0,1]`$, then both $`\mu (x)`$ and $`\sigma (x)`$ are uniquely determined on $`[0,1]`$.

Proof. The proof follows by induction. We define $`y_M=h(x_M)\epsilon =h(\eta )`$. Let us now define $`f_M(y)`$ as in (13) on $`(0,y_{M1})`$ but with boundary condition $`f_M(y_M)=1`$. This implies that $`u_M(y_M;\lambda )=w(x_M;\lambda )`$ where $`u_M(y;\lambda )`$ solves the Schrödinger equation

$$\begin{array}{c}u_M^{\prime \prime }(y)q(y)u_M(y)=\lambda u_M(y),0<y<y_{M1}\hfill \\ u_M(y_{M1})=\frac{1}{f(y_{M1})}>0,\frac{u_M(0;\lambda )}{y}=0.\hfill \end{array}$$
(31)
We have seen in Theorem 3 that $`q(y)`$ can then be uniquely reconstructed on $`(0,y_{M1})`$ since it is known on $`(0,y_M)`$. We can now apply the rest of the proof of Theorem 2 on the interval $`(0,x_{M1})`$. This allows us to reconstruct $`\mu `$ and $`\sigma `$ uniquely on $`(0,x_{M1})`$. Similarly the knowledge of $`\mu `$ and $`\sigma `$ on $`(0,x_m)`$ allows us to reconstruct both coefficients uniquely on $`(0,x_{m1})`$ and the same knowledge on $`(0,x_0)`$ allows reconstruction on $`(0,1)`$. This concludes the proof of the theorem.

## 7 Implementation of the reconstruction

Theorems 1, 2, and 4 provide an explicit method to reconstruct the drift and diffusion coefficients $`\mu `$ and $`\sigma `$ from the measurements. This method is based on three steps. The first step consists of finding the poles of meromorphic functions of the form $`\lambda u(y;\lambda )`$ at fixed $`y`$ from its knowledge for $`\lambda >0`$. In the second step, the potential $`q(y)`$ is reconstructed from the knowledge of the spectral data provided in step 1. Finally $`\mu (x)`$ and $`\sigma (x)`$ are reconstructed from $`q(y)`$ by solving first-order ordinary differential equations.

Step 3 is well-posed in the sense that noise in the potential $`q(y)`$ are at most linearly amplified in the reconstruction of $`\mu (x)`$ and $`\sigma (x)`$. Step 2 is also well posed as the reconstruction of the potential $`q(y)`$ from the spectra is stable, for instance in the sense that an error of order $`\epsilon `$ on the $`l^2`$ norm of $`a_k`$ in (25) implies an error also of order $`\epsilon `$ in the reconstruction of $`q(x)`$ in the $`L^2`$ sense .

Step 1, however, is severely ill-posed in the sense that order-$`\epsilon `$ errors in the measured generate errors greater than order $`\epsilon ^\alpha `$ (for all $`\alpha >0`$) in the $`l^2`$ norm of the coefficients $`a_k`$ in (25). This is a similar result to the reconstruction of diffusion coefficients from boundary measurements, which is known to be severely ill-posed . We can show a Hölder-type stability result in the reconstruction of each eigenvalue. This means that for each eigenvalue $`\lambda _k`$, there exists $`\alpha _k>0`$ such that errors in the data of order $`\epsilon `$ yield an error of order $`\epsilon ^{\alpha _k}`$ on the reconstruction. Unfortunately, the prefactor in front of $`\epsilon ^{\alpha _k}`$ grows exponentially in $`k`$ so that at most only approximately $`|\mathrm{log}\epsilon |`$ eigenvalues can be reconstructed from oder-$`\epsilon `$ accurate measured data. This behavior is a consequence of the results on the meromorphic prolongation of functions with poles, as shown by K. Miller . Since we do not attack the numerical reconstruction here we briefly sketch the result and how it may be used to derive the above stability estimates.

Let us consider the open disk $`D`$ of radius $`R`$ and center $`0`$ in the complex plane and a meromorphic function $`f(\lambda )=u(y_0;\lambda )`$ inside that disk with $`N`$ poles, which we assume lie on the segment $`(R,0)`$. We want to reconstruct these $`N`$ poles from the measurement of the meromorphic function on the segment $`\mathrm{\Gamma }=(0,R)`$. Let us define by $`w(z)`$ the function that takes the value $`0`$ on $`D`$, the value $`1`$ on $`\mathrm{\Gamma }`$, and is harmonic on $`D\backslash \mathrm{\Gamma }`$. One verifies that $`0<w(z)<1`$ on $`D\backslash \mathrm{\Gamma }`$. Let us recast our meromorphic function as $`f(\lambda )=F(\lambda )/_{k=1}^N(\lambda \lambda _k)`$ where $`F(\lambda )`$ is analytic on $`D`$. We assume that

$$|FBh|<\epsilon |B|\text{ on }\mathrm{\Gamma },\text{ and }|F|<E,\text{ on }D.$$

Here $`h(\lambda )`$ are our measurements. The above relation indicates what we mean by error in the measurements and essentially says that $`fh`$ is bounded by $`\epsilon `$ in the supremum norm. Upon slightly modifying the proof of Theorem 1 in , we obtain that it is possible to reconstruct the poles $`\lambda _k`$ of $`f`$ on $`D`$ with an error of order

$$|\delta \lambda _k|\frac{C^NR^N(\epsilon R^N)^{w(\lambda _k)}E^{1w(\lambda _k)}}{|W_k|\left(\underset{lk}{}|\lambda _l\lambda _k|\right)^2}.$$

Here we have bounded $`B`$ by $`R^N`$ on $`D`$. In the case of interest in this paper, we have that

$$RN^2,\lambda _kk^2,\text{and }W_k=u(y^{})\varphi _k(y_0)\varphi _k^{}(y^{}).$$

The radius is chosen as roughly $`1/2(\lambda _N+\lambda _{N+1})`$. At the boundary $`D`$ of the domain, we observe that $`F`$ is bounded by (using Stirling’s approximation)

$$E\underset{lk}{}|\lambda _l\lambda _k|(N!)^2N^{2N+1}e^N.$$

Using these values we find that

$$|\delta \lambda _k|C^N\epsilon ^{w(\lambda _k)}$$

for some positive constant $`C`$. The above constraint is sharp . This shows that the reconstruction of any eigenvalue is Hölder continuous, where the order of continuity $`\alpha =w(\lambda _k)`$. This is therefore a mildly ill-posed problem.

However, although the reconstruction of each eigenvalue is mildly ill-posed, collectively the reconstruction of the whole spectrum is severely ill-posed if one wants an approximation in the $`l^2`$ sense. Indeed for a given noise level $`\epsilon `$, the above formula shows that we can only reconstruct a number of eigenvalues that is proportional to $`|\mathrm{log}\epsilon |`$, for the reconstructed information becomes useless as $`|\delta \lambda _k|1`$ and the corresponding coefficients $`a_k`$ in (25) are not reconstructed at all. Now assuming that $`w(\lambda _k)=1`$, which is the best we can expect, we observe that $`C^N\epsilon =1`$ for $`N|\mathrm{ln}\epsilon |`$. This shows that the reconstruction of the potential $`q(y)`$ is severely ill-posed. Indeed let us assume that the asymptotic coefficient of $`q`$ in (25) are of the form $`a_kk^M`$ for some $`M`$. The $`l^2`$ error on the Fourier coefficients obtained by truncating all coefficients of index greater than $`|\mathrm{ln}\epsilon |`$ is then of order $`N^{M+1/2}|\mathrm{log}\epsilon |^{M+1/2}\epsilon ^\alpha `$ for all $`\alpha >0`$.

## 8 Summary and Conclusions

We have shown that if $`\mu `$ and $`\sigma `$ are known in say, $`x[0,x_b]`$, with $`x_b=h^1(y^{}/2)`$ (Eq. 18), one of them can be reconstructed from a single first-exit time distribution measurement at $`x_0\mathrm{\Lambda }`$. The procedure is formally carried out by transforming the backward Kolmogorov equation for the first-exit time distribution $`w(x;\lambda )`$ into a Schrödinger equation in $`u(y;\lambda )`$ with a potential $`q(y)`$ that is a functional of $`\mu (x)`$ and $`\sigma (x)`$. If this potential is known in the domain $`0yy_bh(x_b)`$, $`q(y)`$ can be reconstructed from the single eigenvalue spectrum obtained by analytically continuing the solution $`u(y;\lambda )`$ (in $`\lambda `$) and finding the poles. It is this pole-finding procedure that is severely ill-posed. The asymptotic limit of the eigenvalues also determine the extent of the transformed domain $`y^{}`$. From the reconstructed $`q(y)`$, one can reconstruct $`\mu (x)`$ or $`\sigma (x)`$, only if the other (or $`\mu (x)/\sigma (x)`$) is also known. Finally, we have shown that the measurement process and the determination of spectra can be repeated a sufficient number of times to obtain $`q(y)`$ (Theorem 3) and $`\mu `$ or $`\sigma `$ (Theorem 4) on the entire domain provided $`q(y)`$ and $`\sigma `$ or $`\mu `$ is initially known on a small interval $`(0,\epsilon )`$.

Our analysis in section 7 shows that the proposed inversion is severely ill-posed even if the hypotheses of the aforementioned theorems are satisfied. However the reconstruction of partial spectral data is mildly ill-posed so we believe the proposed method can be implemented. Let us conclude with the following remark. Our method is based on reconstructing the eigenvalues $`\lambda _k`$ from the knowledge of $`u(x_0;\lambda )`$. The latter may however provide additional information. For instance we observe from (21) that the residuals $`\varphi _k(y)\varphi _k^{}(y^{})u(y^{})`$ can also be reconstructed from the measured data. It is known that the measurement of one spectrum $`\{\lambda _k\}`$ may be sufficient to uniquely reconstruct $`q(y)`$ on the whole interval $`(0,y^{})`$ provided that additional normalization conditions be also measured on the eigenfunctions . It is not clear whether additional information on $`\varphi _k(y)\varphi _k^{}(y^{})u(y^{})`$ may be useful. Preliminary numerical simulations seem to show that reconstructions from a single measurement only allows us to reconstruct part of the potential.

## Acknowledgments

This work was supported by the National Science Foundation (US) through grants DMS-0239097 (GB) and DMS-0206733 (TC) and by the Office of Naval Research through grant N00014-02-1-0089 (GB). GB also acknowledges support from an Alfred P. Sloan Fellowship.

## References


text:
# First-principle density-functional calculation of the Raman spectra of BEDT-TTF

First-principles density-functional calculation of the Raman spectra of BEDT-TTF

K. Brake1, B.J. Powell1, R.H. McKenzie1, M.R. Pederson2 and T. Baruah2,3

1. Department of Physics, University of Queensland, Brisbane, Queensland 4072, Australia

2. Center for Computational Materials Science, Code 6390, Naval Research Laboratory, Washington, D.C. 20375, USA

3. Department of Physics, Georgetown University, Washington, D.C. 20057, USA

Abstract. We present a first-principles density-functional calculation for the Raman spectra of a neutral BEDT-TTF molecule. Our results are in excellent agreement with experimental results. We show that a planar structure is not a stable state of a neutral BEDT-TTF molecule. We consider three possible conformations and discuss their relation to disorder in these systems.

Crystals based on the BEDT-TTF molecule of the form (BEDT-TTF)2X, where X is a monovalent anion such as Cu\[N(CN)2\]Cl, are of great interest both experimentally and theoretically. This is due to the variety of ground states found, and the sensitivity of the ground state to subtle changes in pressure and chemistry . It has been proposed that a minimal model for such crystals is the Hubbard model on an anisotropic triangular lattice , for which we ultimately aim to extract the parameters from quantum chemistry calculations, using NRLMOL, a Gaussian based, all electron density functional theory based electronic structure method developed by Pederson *et al.* .

We have calculated the Raman spectra of a neutral BEDT-TTF molecule. We have assumed the molecule is symmetric about the $`x`$, $`y`$ and $`z`$ axes, that is, that it is in a planar structure, shown in figure 2. The results are shown in figure 1. The range of the plot was chosen to aid comparison with the experimental data from Eldridge *et al.*  reproduced in figure 1. Outside the range of frequencies shown in figure 1 the only notable features are several strong hydrogen stretch peaks around 3000 cm-1 that were also reported by Eldridge *et al.* . There were also several peaks present with imaginary frequencies, which show that the planar structure is not a stable state of the BEDT-TTF molecule. This agrees with the results of quantum chemistry calculations by Demiralp and Goddard .

Figure 1: Left: Raman spectra for BEDT-TTF at 300 K calculated using the generalised gradient approximation form of density functional theory methods. The peaks are broadened by a Gaussian curve with a full width half maximum of 1 cm-1 to aid comparison with experiment. Right: experimental Raman spectra for BEDT-TTF at room temperature from Eldridge *et al.* .

As our calculations have shown that the planar structure is not stable, we have performed calculations on several other conformations of the BEDT-TTF molecule with no symmetry enforced. Two such structures of these are shown in figure 2. The total energy of the two alternative conformations are similar to each other at about 1 eV below that of the planar structure.

It has been shown that conformational changes in the terminal ethylene groups of BEDT-TTF can lead to intrinsically non-magnetic disorder . In both $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Br and $`\beta `$-(BEDT-TTF)2I3 this disorder has been shown to lead to a suppression of the superconducting critical temperature . Understanding this phenomena and the details of the conformational change involved may therefore be important in identifying the pairing symmetry of the superconducting state of these materials . It has been proposed by Ravy *et al.*  that the conformational changes can be modeled as a double harmonic well. Disorder becomes frozen in below T $``$ 80 K in $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Br and $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Cl . It has also been found that the activation energy for terminal ethylene disordering is $``$ 3200 K in $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Br and $``$ 2650 K for $`\kappa `$-(BEDT-TTF)2Cu\[N(CN)2\]Cl . The calculated energy difference between the W and Z structures (see figure 2) is 750 K. Clearly the potential barrier between the relevant conformations then must be greater than the energy difference between the two conformations. However, the activation energy will be determined by the highest energy state which the molecule has to pass through to move from one configuration to another. In this case of a BEDT-TTF molecule moving between the W and Z configurations this may well be close to the the energy difference between the W conformation and the planar conformation, 21200 K, which is an order of magnitude larger than the measured activation energy. On the other hand, it is possible that there is path between the W and Z configurations that does not involve moving through the planar configuration and has a lower activation energy. Also it is clear from the difference between the activation energies of the two salts that interactions with the local environment can change the activation energy of conformational changes. Therefore, although although it seems unlikely that the W and Z conformations are responsible for terminal ethylene disorder in BEDT-TTF salts, this cannot be strictly ruled out.

Figure 2: The conformations of the BEDT-TTF molecule considered in this paper and the relative energies and HOMO-LUMO gaps.

In conclusion, we have calculated the Raman spectra of a neutral BEDT-TTF molecule using first-principles density functional theory methods. Our calculated spectra agrees well with experimental results. The presence of peaks with imaginary frequency show that the structure enforced on the molecule for the calculation is not the stable state. We have shown that the energy difference between the W and Z conformations is not inconsistent with the reported activation energy of $`\kappa `$-(BEDT-TTF)Cu\[N(CN)2\]Br and $`\kappa `$-(BEDT-TTF)Cu\[N(CN)2\]Cl, but a crude estimate of the activation energy between these conformations is an order of magnitude larger than the experimentally observed activation energy.

Acknowledgements

We would like to thank Jens Müller for useful conversations. The work at the University of Queensland was supported by the Australian Research Council. KB was supported in part by US ONR-IFO N00014-03-1-4116. BJP was supported in part by US ONR-IFO N00014-03-1-4115. MRP was supported in part by ONR and the DoD HPC CHSSI initiative. TB was supported in part by US ONR N00014-02-1-1046.