text:
# Response to Compagno and Persico*footnote **footnote *This article is written by V Hnizdo in his private capacity. No official support or endorsement by the Centers for Disease Control and Prevention is intended or should be inferred.

Abstract
Some mathematically incorrect claims of Compagno and Persico in their reply (2002 J. Phys. A: Math. Gen. 35 8965) to my comment on their recent paper on self-dressing and radiation reaction in classical electrodynamics are pointed out.

Compagno and Persico (CP) have replied  to my comment  on their paper on self-dressing and radiation reaction in classical electrodynamics . CP acknowledge the main point of the comment, namely that the expression for the time-averaged electromagnetic self-force obtained in  for the test charge of a Bohr–Rosenfeld field-measurement procedure and rejected in  as incorrect can be obtained also using a formula for the self-force which they derived by different means in . In view of this fact, CP now endorse the expression in question as correct. However, some claims in their reply call for my response.

The expression for the time-averaged self-force on a spherical uniform charge $`q`$ of radius $`a`$ obtained in  reads

$$\overline{F}=\frac{q^2}{T}_0^Tdt^{}Q(t^{})f(t^{})$$
(1)
where

$$f(t^{})=\frac{1}{2a^3}(2\chi )(22\chi \chi ^2)\mathrm{\Theta }(2\chi )\chi =\frac{Tt^{}}{a}.$$
(2)
Here, the speed of light $`c=1`$ and $`Q(t)`$ is the charge’s one-dimensional trajectory, which is subject to the conditions that $`Q(t)=0`$ for $`t0`$ and $`|Q(t)|a`$, $`|\mathrm{d}Q(t)/\mathrm{d}t|c`$ for $`0<t<T`$; $`\mathrm{\Theta }(x)`$ is the Heaviside step function. Instead of the simple closed-form expression (2) for the function $`f(t^{})`$, CP counter-proposed in  the expression (normalized here to conform with (2)):

$$f(t^{})=\frac{2}{3V^2}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!}[\delta ^{(n+1)}(Tt^{})\delta ^{(n+1)}(t^{})]r^{n1}V=\frac{4}{3}\pi a^3$$
(3)
where

$$r^{n1}=_{|𝒓_1|<a}d𝒓_1_{|𝒓_2|<a}d𝒓_2|𝒓_2𝒓_1|^{n1}=\frac{72V^2(2a)^{n1}}{(n+5)(n+3)(n+2)}.$$
(4)
(The above closed-form expression for $`r^{n1}`$ was given subsequently in .) My first point here has to be unfortunately of a rather trivial nature. The sign of expression (3) is correct, if, as in the standard notation of elementary calculus, $`u^{(n)}[v(x)]\mathrm{d}^nu(y)/\mathrm{d}y^n|_{y=v(x)}`$. Thus, e.g., $`\delta ^{(n)}(Tt)\mathrm{d}^n\delta (x)/\mathrm{d}x^n|_{x=Tt}`$. With this notation, there is no sign misprint in the function $`f(t^{})`$ in  that CP now want to correct—but if $`\delta ^{(n)}(Tt)`$ meant instead $`\mathrm{d}^n\delta (Tt)/\mathrm{d}t^n=(1)^n\mathrm{d}^n\delta (x)/\mathrm{d}x^n|_{x=Tt}`$, as CP suggest, then not only the sign but also the factor $`(1)^n`$ would be there in error. In either case, the function $`f(t^{})`$ is now given incorrectly by the expression (3) of —it either has the wrong overall sign, or the factor $`(1)^n`$ there must be omitted.

Contrary to an assertion of CP, I have never claimed that the expression (3) for the function $`f(t^{})`$ is incorrect. I have rather pointed out in  that, in order to obtain the time-averaged self-force (1), this expression was used incorrectly by CP in the requisite integration. Since this is an integration with finite limits involving high-order derivatives of the delta function, it cannot be performed as simplistically as CP have done, obtaining

$$\overline{F}_{\mathrm{CP}}=\frac{2\rho ^2}{3T}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!}[Q^{(n+1)}(T)Q^{(n+1)}(0)]r^{n1}\rho =\frac{q}{V}.$$
(5)
This is the ‘exact’ expression for the time-averaged self-force $`\overline{F}`$ that CP still claim in  to be correct, despite the fact that it is inconsistent with their newly adopted approval of the expression (2) for the function $`f(t^{})`$. This can be shown by evaluating the time-averaging integral (1) for $`\overline{F}`$ using the expression (2) for $`f(t^{})`$ and the Taylor expansion of the trajectory $`Q(t)`$ about the point $`t=T`$:

$`\overline{F}`$ $`=`$ $`{\displaystyle \frac{q^2}{T}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{n!}}Q^{(n)}(T){\displaystyle _0^T}dt^{}(Tt^{})^nf(t^{})`$ (6)
$`=`$ $`{\displaystyle \frac{q^2}{T}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^n}{n!}}Q^{(n)}(T)\left\{\alpha _n\left[a^{n2}\left({\displaystyle \frac{2\kappa ^{n+1}}{n+1}}{\displaystyle \frac{3\kappa ^{n+2}}{n+2}}+{\displaystyle \frac{\frac{1}{2}\kappa ^{n+4}}{n+4}}\right)+\alpha _n\right]\mathrm{\Theta }(2\kappa )\right\}`$ (7)
where

$$\alpha _n=\frac{48(2a)^{n2}n}{(n+1)(n+2)(n+4)}\kappa =\frac{T}{a}.$$
(8)
This is the correct expression for the time-averaged self-force in terms of the derivatives $`Q^{(n)}(T)lim_{tT^{}}\mathrm{d}^nQ(t)/\mathrm{d}t^n`$. (Note for completeness that the correct expression for $`\overline{F}`$ in terms of the derivatives $`Q^{(n)}(0)lim_{t0^+}\mathrm{d}^nQ(t)/\mathrm{d}t^n`$ is given by equation (25) of .) As $`\alpha _n=\frac{2}{3}nr^{n2}/V^2`$, using (7) we can write $`\overline{F}`$ for $`T2a`$ as

$$\overline{F}=\frac{q^2}{T}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!}Q^{(n)}(T)\alpha _n=\frac{2\rho ^2}{3T}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n}{n!}Q^{(n+1)}(T)r^{n1}T2a$$
(9)
which happens to be equal to the first of the two terms in the expression (5) of CP for $`\overline{F}`$. This demonstrates that expression (5) is incorrect; for $`T2a`$, (5) can be corrected by dropping its second term, but for $`T<2a`$, no such simple correction is possible.

Far from being ‘convenient’ for an ‘exact’ evaluation of the time-averaged self-force $`\overline{F}`$, formula (3) is a purely formal expression that has no practical application in an integration with finite limits. Its use by CP has led to the erroneous expression (5) for $`\overline{F}`$; when CP use it in  to prove the equivalence of expressions (2) and (3), they revert the Taylor expansions $`_n(1)^n\delta ^{(n+1)}(t)r^n/n!`$ to $`\delta ^{}(tr)`$ before performing the finite-limit integration.

So far, CP have responded to my criticism  of their re-analysis  of the Bohr–Rosenfeld field-measurement procedure only by making mathematically incorrect claims. All these involved rather simple mathematical points about which there should have been no need of explicating since my paper  of 2000 (where I derived expression (2) using elementary calculus of the delta function), if not already since my comment  of 1999 (where I used Fourier transform methods). It is regrettable that such points have deflected from the interesting issues of physics relating to the famous Bohr–Rosenfeld analysis.


format: text

text:
# Intelligence and Cooperative Search by Coupled Local Minimizers

## 1 Introduction

A large variety of problems arising in engineering, physics, chemistry and economics can be formulated as optimization problems, either constrained or unconstrained having continuous or discrete variables. A well developed area in optimization theory are local optimization methods including conjugate gradient, quasi-Newton, Levenberg-Marquardt, sequential quadratic programming etc. \[Bertsekas, 1996; Fletcher, 1987; Gill et al., 1981\]. Continuous-time optimization methods have been developed within the area of neural networks \[Cichocki & Unbehauen, 1994\]. Popular methods for global exploration of the search space are e.g. simulated annealing \[Kirkpatrick et al., 1983\] and genetic algorithms \[Goldberg, 1989\]. Within standard local optimization techniques one often applies multi-start local optimization, by trying different starting points and running the processes independently from each other and selecting the best result from all trials. On the other hand, for training of neural networks such as multilayer perceptrons (MLP) it is well-known that instead of training several MLPs for random choices of small initial weights and selecting the best of all trained networks, one better forms a committee network which is based upon all trained networks \[Arbib, 1995; Bishop, 1995\]. In other words, the training efforts of the less optimal networks are not entirely useless but can be employed in order to improve the estimate in view of the bias-variance trade-off \[Bishop, 1995\]. Unfortunately, committee networks are only applicable to static nonlinear regression and classification problems.

In this paper we propose a new methodology of coupled local minimizers (CLM) for solving continuous nonlinear optimization problems. We pose a somewhat similar challenge as for committee networks but within a different and broader context of solving differentiable optimization problems. The aim is to (on-line) combine the results from local optimizers in order to let the ensemble generate a local minimum that is better than the best result obtained from all individual local minimizers. We show how improved local minima can be obtained by having interaction and information exchange between the local search processes. This is realized through state synchronization constraints that are imposed between the local minimizers by incorporating principles of master-slave dynamics. Synchronization theory has been intensively studied within the area of chaotic systems and secure communications \[Chen & Dong, 1998; Pecora & Carroll, 1990; Suykens et al., 1996, 1997, 1998; Wu & Chua, 1994\]. The CLM method is related to Lagrange programming network approaches for chaos synchronization \[Suykens & Vandewalle, 2000\], where identical or generalized synchronization constraints are imposed on dynamical systems. CLMs also fit within the framework of Cellular Neural Networks (CNN) \[Chua & Roska, 1993; Chua et al., 1995; Chua, 1998\]. By considering the objective of minimizing the average cost of an ensemble of local minimizers subject to pairwise synchronization constraints, a continuous-time optimization algorithm is studied according to Lagrange programming networks \[Cichocki & Unbehauen, 1994; Zhang & Constantinides, 1992\]. The resulting continuous-time optimization algorithm is described by an array of coupled nonlinear cells or a one-dimensional CNN with bi-directional coupling.

We show how to obtain intelligence from CLMs. This is done by considering a problem of static nonlinear regression with MLPs from given noisy measurement data. In this case CLMs correspond to coupled backpropagation processes \[Rumelhart et al., 1986; Werbos, 1990\]. In order to obtain an MLP with good generalization performance (i.e. an intelligent solution) it is usually needed to consider a regularization term in addition to the original sum squared error cost function (fitting error) which is defined on the training data \[Bishop, 1995; MacKay, 1992; Poggio & Shelton C. 1999; Suykens & Vandewalle, 1998\]. When applying CLMs such a regularization term is not needed. It turns out that the initial distribution of the initial states of the several local minimizers can play this role, such that one only has to optimize the training set error. This is consistent with insights from Bayesian learning of neural networks where the regularization term is related to the choice of the prior on the unknown interconnection weights, which expresses that the initial weights should be chosen small (weight decay).

The role of the initial CLM state is not only important for neural networks but also with respect to potential energy surface optimization of Lennard-Jones clusters, a second class of problems that is investigated in this paper. It is considered to be an important benchmark problem in the area of protein folding \[Šali et al., 1994; Neumaier, 1997; Wales & Scheraga, 1999\]. We show that CLMs are able to detect the global minimum of (LJ)38 which possesses a double-funnel energy landscape and is known to be a challenging test-case \[Wales & Doye, 1997; Wales et al., 1998; Wales & Scheraga, 1999\]. The CLM method has also been successfully applied to larger scale problems. The cooperative search mechanism in CLMs is obtained by solving a linear programming problem in the unknown soft constraint penalty factors. In this way a maximal decrease in the average energy cost of the ensemble is achieved.

This paper is organized as follows. In Section 2 we introduce coupled local minimizers. In Section 3 we discuss how to obtain optimal interaction between local minimizers. In Section 4 we present new insights on coupled backpropagation processes and intelligence. In Section 5 we apply coupled local minimizers to the optimization of Lennard-Jones clusters.

## 2 A CNN of Coupled Local Minimizers

Consider the minimization of a twice continuously differentiable cost function $`U(𝐱)`$ with $`𝐱^n`$. Let us take an ensemble consisting of $`q`$ local minimizers. We aim at minimizing the average energy cost $`U=\frac{1}{q}_{i=1}^qU[𝐱^{(i)}]`$ of this ensemble subject to pairwise state synchronization of the particles:

$$\underset{𝐱^{(i)}^n}{\mathrm{min}}U\mathrm{such}\mathrm{that}𝐱^{(i)}𝐱^{(i+1)}=0,i=1,2,\mathrm{},q$$
(1)
with particle states $`𝐱^{(i)}^n`$ for $`i=1,2,\mathrm{},q`$ and boundary conditions $`𝐱^{(0)}=𝐱^{(q)}`$, $`𝐱^{(q+1)}=𝐱^{(1)}`$. The synchronization constraints have to be achieved in an asymptotical sense, i.e. the particles have to reach the same final state. One defines the augmented Lagrangian

$$(𝐱^{(i)},\lambda ^{(i)})=\frac{\eta }{q}\underset{i=1}{\overset{q}{}}U[𝐱^{(i)}]+\frac{1}{2}\underset{i=1}{\overset{q}{}}\gamma _i𝐱^{(i)}𝐱^{(i+1)}_2^2+\underset{i=1}{\overset{q}{}}\lambda ^{(i)},[𝐱^{(i)}𝐱^{(i+1)}]$$
(2)
with Lagrange multipliers $`\lambda ^{(i)}^n`$ $`(i=1,2,\mathrm{},q)`$. The different terms in this Lagrangian are the objective function, the soft constraints and the hard constraint related to the pairwise state synchronization constraints. The penalty factors $`\gamma _i`$ emphasize the importance of each of the soft synchronization constraints. From this augmented Lagrangian one obtains the Lagrange programming network \[Zhang & Constantinides, 1992\]

$$\{\begin{array}{ccc}\hfill \dot{𝐱}^{(i)}& =& _{𝐱^{(i)}}(𝐱^{(i)},\lambda ^{(i)})\hfill \\ \hfill \dot{\lambda }^{(i)}& =& _{\lambda ^{(i)}}(𝐱^{(i)},\lambda ^{(i)}),i=1,2,\mathrm{},q.\hfill \end{array}$$
(3)
This is basically a continuous-time optimization algorithm for solving the given constrained optimization problem. One obtains the following array of coupled local minimizers (CLM)

$$\{\begin{array}{ccc}\hfill \dot{𝐱}^{(i)}& =& \frac{\eta }{q}_{𝐱^{(i)}}U[𝐱^{(i)}]+\gamma _{i1}[𝐱^{(i1)}𝐱^{(i)}]\gamma _i[𝐱^{(i)}𝐱^{(i+1)}]+\lambda ^{(i1)}\lambda ^{(i)}\hfill \\ \hfill \dot{\lambda }^{(i)}& =& 𝐱^{(i)}𝐱^{(i+1)},i=1,2,\mathrm{},q\hfill \end{array}$$
(4)
with learning rate $`\eta `$. This array consists of a number of $`q`$ coupled nonlinear cells and is considered to be a special case of cellular neural networks (CNN) \[Chua, 1998\].

The basic mechanism of state synchronization with coupled local minimizers is illustrated in Fig. 1 for a double potential well problem with two particles that are searching for the global minimum. The main idea here is to impose that they should reach the same final position. A simplified CLM considered in this case is

$$\{\begin{array}{ccc}\hfill \dot{x}& =& _xU(x)(xz)\lambda \hfill \\ \hfill \dot{z}& =& _zU(z)+(xz)+\lambda \hfill \\ \hfill \dot{\lambda }& =& xz\hfill \end{array}$$
(5)
where $`x,z,\lambda `$. This is derived from the augmented Lagrangian

$$(x,z,\lambda )=\frac{\eta }{2}[U(x)+U(z)]+\frac{1}{2}(xz)^2+\lambda (xz).$$
(6)
As a result of the constraining of the system, the two particles are enforced to take a decision about which valley to choose. When the initial states of the two particles are located in different valleys, one observes that they are always reaching the global minimum. This phenomenon is independent from the steepness and shape of the valleys. In this process, the search space is in fact duplicated and the local optimizers are exchanging information through the synchronization constraint.

Note that it is not necessarily guaranteed beforehand that the equality constraints in (1) will be exactly realized by the CLM, as is known for Lagrange programming networks. In this case one could consider as CLM goal that the state synchronization constraints should be realized in an approximate sense, i.e. bringing the states close to each other. Eventually, the equality constraints could be replaced by a set of inequality constraints which would implement this.

## 3 CLMs and Optimal Cooperative Search

In this Section, we discuss how to obtain an optimal interaction between the several local minimizers. This is done by making a suitable choice of the penalty weights $`\gamma _i`$. In this design procedure we aim at maximally decreasing the instantaneous average energy cost of the ensemble

$$\begin{array}{ccc}\hfill \frac{dU}{dt}& =& \frac{1}{q}\underset{i=1}{\overset{q}{}}\frac{U[𝐱^{(i)}]}{𝐱^{(i)}},\dot{𝐱}^{(i)}\hfill \\ & & \\ & =& \frac{1}{q}\underset{i=1}{\overset{q}{}}\frac{U[𝐱^{(i)}]}{𝐱^{(i)}},\frac{\eta }{q}\frac{U[𝐱^{(i)}]}{𝐱^{(i)}}+\gamma _{i1}[𝐱^{(i1)}𝐱^{(i)}]\hfill \\ & & \gamma _i[𝐱^{(i)}𝐱^{(i+1)}]+\lambda ^{(i1)}\lambda ^{(i)}.\hfill \end{array}$$
(7)
For a given state vector $`𝐱=[𝐱^{(1)};𝐱^{(2)};\mathrm{};𝐱^{(q)}]`$ and costate of the ensemble $`\lambda =[\lambda ^{(1)};\lambda ^{(2)};\mathrm{};`$ $`\lambda ^{(q)}]`$ the latter expression is affine in $`\gamma =[\gamma _1;\gamma _2;\mathrm{};\gamma _q]`$ and the interaction between the local minimizers is optimized then by solving the linear program (LP):

$$\underset{\gamma ^q}{\mathrm{min}}\frac{dU}{dt}_{𝐱,\lambda }\mathrm{such}\mathrm{that}\underset{¯}{\gamma }<\gamma _i<\overline{\gamma },i=1,2,\mathrm{},q$$
(8)
where $`\underset{¯}{\gamma },\overline{\gamma }`$ are user-defined lower- and upperbounds. The resulting $`\gamma _i`$ values are kept constant for simulation of the CLM over a certain time interval $`\mathrm{\Delta }T`$. The difference between the $`\gamma _i`$ values in (4) causes a similar effect as is known in master-slave (drive-response) synchronization which has been studied for secure communications using chaos \[Chen & Dong, 1998; Pecora & Carroll, 1990; Suykens et al., 1997, 1998; Wu & Chua, 1994\]. In this way, successfully performing local optimizers impose their state to the other part of the ensemble, which leads to cooperative behaviour.

The step size $`\eta `$ is determined here by imposing a target evolution law

$$d(UU^{})/dt=\alpha (UU^{})$$
(9)
where $`\alpha `$ is a user-defined positive real constant and $`U^{}`$ is an estimate of the energy cost value at the global minimum (alternatively one adds a constant value to the cost function such that the energy cost is guaranteed to be positive everywhere). As a result the step size $`\eta `$ becomes $`\eta =[q_{i=1}^q\frac{U[𝐱^{(i)}]}{𝐱^{(i)}},𝐡+q\alpha (_iU[𝐱^{(i)}]qU^{})]/_{i=1}^q\frac{U[𝐱^{(i)}]}{𝐱^{(i)}},\frac{U[𝐱^{(i)}]}{𝐱^{(i)}}`$ subject to additional user-defined constraints $`\underset{¯}{\eta }<\eta <\overline{\eta }`$ where $`𝐡=\gamma _{i1}[𝐱^{(i1)}𝐱^{(i)}]\gamma _i[𝐱^{(i)}𝐱^{(i+1)}]+\lambda ^{(i1)}\lambda ^{(i)}`$ by definition. The values of $`\gamma _i,\eta `$ are re-scheduled each time in between simulations of the array over user-defined time intervals $`\mathrm{\Delta }T`$.

According to the idea of small-world networks \[Watts & Strogatz, 1998\], the information spreading within the ensemble can be further improved, e.g. by re-numbering a part of the ensemble (state vectors together with their corresponding co-state) at random every $`c\mathrm{\Delta }T`$ steps where $`c_0`$. In this way the pairwise synchronization constraints vary among the different local minimizers during optimization.

In order to obtain insight in the influence of the CLM tuning parameters let us consider the optimization of the cost function

$$\begin{array}{c}U(𝐱)=\frac{1}{2n}a\underset{i=1}{\overset{n}{}}x_i^2+8n4n\underset{i=1}{\overset{n}{}}\mathrm{cos}(\omega _1x_i)4n\underset{i=1}{\overset{n}{}}\mathrm{cos}(\omega _2x_i)\end{array}$$
(10)
with $`a=0.01`$, $`\omega _1=0.2`$, $`\omega _2=1`$ in a $`n=10`$ dimensional search space \[Styblinski & Tang, 1990\]. The global minimum for this problem is known and located at the origin with $`U(0)=0`$ and $`U^{}=0`$ is taken. The solutions from CLMs improve with respect to multi-start local optimization, even when a quasi-Newton method instead of steepest descent optimizers (as for the CLM) are taken. In these comparisons the initial states have been random uniformly generated in $`[\mathrm{20\hspace{0.17em}\hspace{0.17em}20}]^n`$ where the same initial states have been taken for the different algorithms and several runs of the algorithms have been done. CLMs with $`q20`$ are reaching the global minimum in this case. From this example it follows that in order to optimize a surface over a certain region in search space, the number $`q`$ needed to achieve a good performance depends on the complex shape of the surface or typically on the number of local minima per volume in search space that one intends to explore. Experiments suggest that a factor $`\overline{\gamma }/\underset{¯}{\gamma }=10`$ is usually a good choice. In this way the energy level or cost function of the optimizers remains bundled. Otherwise the simulations could lead to excessive exploration and decreased speed of convergence. Typically, larger values of $`\overline{\gamma }`$ will require shorter $`\mathrm{\Delta }T`$ intervals. The choice of the pressure coefficient $`\alpha `$ will control the average energy decrease of the ensemble. Demanding a high pressure $`\alpha `$, however, might lead to less exploration and worse local minima. The bounds on $`\eta `$ are taken here as $`\underset{¯}{\eta }=10^2`$, $`\overline{\eta }=10^3`$. Although it is not necessary, the CLM performance can be further enhanced by applying the small-world network concept, by re-numbering at random (here 20% of the local minimizers) every $`c\mathrm{\Delta }T`$ steps where $`c=5`$ typically, leading to improved information spreading among the optimizers. For the simulations here and in the sequel a Runge-Kutta integration rule with adaptive step size (absolute and relative error equal to $`10^2`$) has been used.

## 4 Intelligence and Coupled Backpropagation

In this Section we apply CLMs to the training of multilayer perceptron networks

$$y=w^T\mathrm{tanh}(Vu+\beta )$$
(11)
with input vector $`u^m`$, output $`y`$, interconnection weights $`w^{n_h}`$, $`V^{n_h\times m}`$ and bias vector $`\beta ^{n_h}`$, where $`n_h`$ denotes the number of hidden units. Let us denote the unknown parameter vector as $`\theta =[w;V(:);\beta ]`$ which stacks all the weights into one single vector. Since the introduction of the backpropagation algorithm for training of neural networks \[Rumelhart et al., 1986; Werbos, 1990\], important insights have been obtained about regularization (weight decay) of the cost function in order to have a good generalization ability in view of the bias-variance trade-off \[Bishop, 1995; MacKay, 1992; Suykens & Vandewalle, 1998\]. It is well-known that instead of optimizing the cost function

$$\underset{\theta }{\mathrm{min}}J(\theta )=\frac{1}{2}\underset{k=1}{\overset{N}{}}[d_ky_k(\theta )]^2$$
(12)
where $`d_k`$ are the desired target outputs (with given data set $`\{u_k,d_k\}_{k=1}^N`$ of $`N`$ training data), one better considers the following cost function with regularization

$$\underset{\theta }{\mathrm{min}}J^{reg}(\theta ;\mu ,\zeta )=\zeta \frac{1}{2}\underset{k=1}{\overset{N}{}}[d_ky_k(\theta )]^2+\mu \frac{1}{2}\theta ^T\theta $$
(13)
where $`\mu ,\zeta `$ are additional hyperparameters which can be automatically tuned within the framework of Bayesian learning. The second term is related to a prior distribution on the weights which typically expresses that the initial weights are chosen small. The first term is associated with the likelihood. When networks are overparametrized (too many hidden units), regularization enables to work implicitly with a number of effective parameters that is less than the actual number of interconnection weights such that overfitting can be avoided.

Let us now consider a CLM by taking as state vectors $`𝐱^{(i)}`$ ($`i=1,\mathrm{},q`$) the unknown parameter vectors $`\theta ^{(i)}`$. Hence, the CLM consists in fact of coupled backpropagation processes with weight vector synchronization. In Fig. 2 an example of estimating a sine function from given noisy data is given where an overparametrized MLP is taken. Simulation results show that by CLMs one is able to obtain a similar smooth solution as from (13) with Bayesian learning but by optimizing the original cost function (12) without regularization. In order to obtain this intelligent solution (i.e. a good generalization performance), the initial distribution of $`𝐱^{(i)}(0)`$ ($`i=1,\mathrm{},q`$) (at time 0) plays an important role, similar with the choice of $`\mu `$ in (13). Consider the following distribution of the initial CLM state

$$p(𝐱(0))\mathrm{exp}[\frac{1}{2\sigma ^2}𝐱(0)^T𝐱(0)]$$
(14)
which assumes that the states of the local minimizers $`𝐱^{(i)}(0)`$ are normally distributed. The value of $`\sigma `$ then plays the role of a bifurcation hyperparameter. For $`\sigma `$ chosen very large a solution with bad generalization is obtained, while for $`\sigma `$ small (small initial weights) a good generalization performance is obtained. Hence, CLMs can avoid the overfitting phenomenon without the use of a regularization in the sum squared error cost function. In order to obtain the results with good generalization in Fig.2, the following CLM tuning parameters were taken: $`q=20`$, $`\sigma =0.1`$, $`\alpha =10^5`$, $`\underset{¯}{\gamma }=10^6`$, $`\overline{\gamma }=10^7`$, $`\underset{¯}{\eta }=10^2`$, $`\overline{\eta }=10^{10}`$, $`U^{}=0`$ and 20% of the local minimizers is re-numbered after $`c=5`$ with $`\mathrm{\Delta }T=10^9`$.

## 5 Optimization of Lennard-Jones Clusters

In this Section we discuss the application of CLMs to the optimization of Lennard-Jones clusters. In predicting the three-dimensional structure of proteins from amino acid sequences, potential energy surface (PES) minimization is often related to the native structure of the protein. Optimization of Lennard-Jones (LJ) clusters is considered to be an important benchmark problem at this point \[Šali et al., 1994; Neumaier, 1997; Wales & Scheraga, 1999\]. Recently, detailed studies have been reported about optimization of (LJ)N with $`N150`$ by means of basin-hopping techniques in comparison with many other approaches \[Wales & Doye, 1997; Wales et al., 1998; Wales & Scheraga, 1999\]. The cost function in this case is given by

$$U_{LJ}=4\underset{i<j}{}(\frac{1}{r_{ij}^{12}}\frac{1}{r_{ij}^6})$$
(15)
where $`r_{ij}`$ the Euclidean distance between atom $`i`$ and $`j`$ ($`j=1,\mathrm{},N`$). (LJ)38 which possesses a double-funnel energy landscape and is known to be an interesting test-case.

Fig. 3 shows the result of a CLM run applied to (LJ)38. The CLM tuning parameters in this case are $`q=50`$, $`\alpha =10^5`$, $`\underset{¯}{\gamma }=10^6`$, $`\overline{\gamma }=10^7`$, $`\underset{¯}{\eta }=10^2`$, $`\overline{\eta }=10^6`$, $`U^{}=0`$ ($`U_{LJ}^\delta =U_{LJ}+\delta `$ with $`\delta =200`$ is optimized in order to ensure a positive cost function everywhere) and 10% of the local minimizers is re-numbered after $`c=5`$ with $`\mathrm{\Delta }T=10^8`$. The CLM energy bundle caused by state synchronization of the local minimizers is clearly visible. Like in the application of CLMs to the training of neural networks, the choice of $`\sigma `$ in $`p(𝐱(0))`$ plays an important role. In Fig. 3 the differences are illustrated for $`\sigma =0.1`$ versus $`\sigma =1`$. The choice of this prior implicitly corresponds in fact to an additional confining potential energy term, which has been explicitly considered in effective potential minimization methods \[Schelstraete & Verschelde, 1997\] (in a different form). However, the effective potential minimization method has problems with detecting e.g. the global minimum of (LJ)8,9, which on the other hand are easily found by CLMs. In the example of (LJ)38 a classical local optimization technique (quasi-Newton and scaled conjugate gradient for large scale problems) has been used for post-processing and fine-tuning (about 10-20 additional iterations) of the CLM results. The variety of generated solutions includes the global minimum -173.9284.

In the application of CLMs to larger clusters, first a shifted version of the problem $`U_{LJ}^{shift}=4_{i<j}[\frac{1}{(r_{ij}+\mu )^{2\nu }}\frac{1}{(r_{ij}+\mu )^\nu }]`$ has been minimized in order to avoid ill-conditioning in the LP problem and excessively large values in the initial CLM evolution. For (LJ)150 $`\mu =0.1`$, $`\nu =3`$ has been chosen. These results have been taken then as initialization for optimization of $`U_{LJ}`$ and are, eventually, post-processed by a classical local optimization method. The CLM tuning parameters in the (LJ)150 case are $`q=100`$, $`\sigma =0.1`$, $`\alpha =10^5`$, $`U^{}=0`$, $`U_{LJ}^{shift,\delta }=U_{LJ}^{shift}+\delta `$ with $`\delta =3000`$ and $`U_{LJ}^\delta =U_{LJ}+\delta `$ with $`\delta =1000`$. The other tuning parameters for (LJ)150 are the same as for the (LJ)38 example. This CLM run yields a PES value of -872.91, while the best known minimum is -893.31 in this case, obtained by counting nearest-neighbour interactions for icosahedral packing schemes (Northby). However, like simulated annealing and basin-hopping approaches, our focus is here on unbiased methods which are generally applicable.

One can argue that CLMs are providing good solutions given the user-defined initial distribution (prior). In this sense, CLM results can potentially be further improved by exploiting more informative prior knowledge about the problem. The simulation results show that this becomes more important for larger clusters. Simulations have been done on a Sun Ultra 2 workstation in Matlab with cmex implementation of the PES evaluation and gradient calculations. CLM running times are roughly comparable with $`q`$ times applying a standard gradient based local optimization methods (conjugate gradient).

## 6 Conclusions

We have introduced a new optimization method of coupled local minimizers. This is done by considering the average energy cost of the total ensemble subject to pairwise synchronization constraints that are asymptotically reached. The resulting continuous-time optimization algorithm from Lagrange programming networks is an array of coupled nonlinear cells or CNN. For problems of MLP training in supervised learning, we have shown that CLMs are able to produce intelligent solutions, in the sense that one can obtain a good generalization ability without taking a regularization term in the cost function. The obtained insights are consistent with Bayesian learning methods for MLPs. We have also shown how one can let the ensemble of local minimizers cooperate in an optimal way. This is done by solving additional linear programming problems. The CLM method has been successfully applied to the optimization of Lennard-Jones clusters which is an important benchmark problem in the area of protein folding. We expect that the proposed CLM methodology may offer new insights for the application of continuous-time optimization algorithms to NP complete problems in general. It also emphasizes the importance of the role played by the initial state distribution of the local minimizers that form the CLM ensemble.

Acknowledgements. This research work was carried out at the ESAT laboratory and the Interdisciplinary Center of Neural Networks ICNN of the Katholieke Universiteit Leuven, in the framework of the FWO projects Learning and Optimization: an Interdisciplinary Approach and Collective Behaviour and Optimization: an Interdisciplinary Approach, the Belgian Program on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister’s Office for Science, Technology and Culture (IUAP P4-02 & IUAP P4-24) and the Concerted Action Project MEFISTO of the Flemish Community. Johan Suykens is a postdoctoral researcher with the National Fund for Scientific Research FWO - Flanders.

References

Arbib, M.A. (Ed.)  The Handbook of Brain Theory and Neural Networks (MIT Press, Cambridge MA).

Bertsekas, D.P.  Constrained Optimization and Lagrange Multiplier Methods (Athena Scientific).

Bishop, C.M.  Neural networks for pattern recognition (Oxford University Press).

Chen, G. & Dong, X.  From Chaos to Order - Perspectives, Methodologies, and Applications (World Scientific Pub. Co., Singapore).

Chua, L.O. & Roska, T.  “The CNN Paradigm,” IEEE Trans. Circuits and Systems-I, 40(3), 147-156.

Chua, L.O., Hasler, M., Moschytz, G.S. & Neirynck, J.  “Autonomous Cellular Neural Networks: a Unified Paradigm for Pattern Formation and Active Wave Propagation,” IEEE Trans. Circuits and Systems-I, 42(10), 559-577.

Chua, L.O.  CNN: a Paradigm for Complexity (World Scientific Pub. Co., Signapore).

Cichocki, A. & Unbehauen, R.  Neural Networks for Optimization and Signal Processing (Wiley, Chichester).

Fletcher, R.  Practical Methods of Optimization (Wiley, Chichester).

Gill, P.E., Murray, W. & Wright, M.H.  Practical Optimization (Academic Press, London).

Goldberg, D.E.  Genetic Algorithms in Search, Optimization and Machine Learning (Addison-Wesley, Reading Mass.).

Kirkpatrick, S., Gelatt, C.D. & Vecchi, M.  “Optimization by simulated annealing,” Science 220, 621-680.

MacKay, D.J.C.  “Bayesian Interpolation,” Neural Computation 4(3), 415-447.

Neumaier, A.  “Molecular modeling of proteins and mathematical prediction of protein structure,” SIAM Rev. 39, 407-460.

Pecora, L.M. & Carroll, T.L.  “Synchronization in Chaotic Systems,” Phys. Rev. Lett. 64, 821-824.

Poggio, T. & Shelton, C.  “Machine Learning, Machine Vision and the Brain,” AI Magazine 20(3), 37-55.

Rumelhart, D.E., Hinton, G.E. & Williams, R.J.  “Learning representations by back-propagating errors,” Nature 323, 533-536.

Šali, A., Shakhnovich, E. & Karplus, M.  “How does a protein fold?,” Nature 369, 248-251.

Schelstraete, S. & Verschelde, H.  “Finding minimum-energy configurations of Lennard-Jones clusters using an effective potential,” J. Phys. Chem. A 101, 310-315.

Styblinski, M.A. & Tang, T.-S.  “Experiments in nonconvex optimization: stochastic approximation with function smoothing and simulated annealing,” Neural Networks 3, 467-483.

Suykens, J.A.K., Vandewalle, J.P.L. & De Moor, B.L.R.  Artificial Neural Networks for Modelling and Control of Non-Linear systems, (Kluwer Academic Publishers, Boston).

Suykens, J.A.K., Curran, P.F., Vandewalle, J. & Chua, L.O.  “Robust nonlinear H synchronization of chaotic Lur’e systems,” IEEE Trans. Circuits and Systems-I 44, 891-904.

Suykens, J.A.K. & Vandewalle, J. (Eds.)  Nonlinear Modeling: advanced black-box techniques (Kluwer Academic Publishers, Boston).

Suykens, J.A.K., Yang, T. & Chua, L.O.  “Impulsive synchronization of chaotic Lur’e systems by measurement feedback,” Int. J. Bifurcation and Chaos 8, 1371-1381.

Suykens, J.A.K. & Vandewalle, J.  “Chaos synchronization: a Lagrange programming network approach,” Int. J. Bifurcations and Chaos (special issue on chaos control and synchronization), 10(4), 797-810.

Wales, D.J. & Doye, J.P.K.  “Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms,” J. Phys. Chem. A 101, 5111-5116.

Wales, D.J., Miller, M.A. & Walsh, T.R.  “Archetypical energy landscapes,” Nature 394, 758-760.

Wales, D.J. & Scheraga, H.A.  “Global optimization of clusters, crystals, and biomolecules,” Science 285, 1368-1372.

Watts, D. J. & Strogatz, S.H.  “Collective dynamics of ’small-world’ networks,” Nature 393, 440-442.

Wu, C.W. & Chua, L.O.  “A unified framework for synchronization and control of dynamical systems,” Int. J. Bifurcation and Chaos 4, 979-989.

Werbos, P.  “Backpropagation through time: what it does and how to do it,” Proceedings of the IEEE 78(10), 1150-1560.

Zhang, S. & Constantinides, A.G.,  “Lagrange programming neural networks,” IEEE Trans. Circuits and Systems-II 39, 441-452.

Captions of Figures

Fig. 1 Illustration of the basic state synchronization mechanism in coupled local minimizers (CLM): (A) double well cost function $`U(x)=x^416x^2+5x+100`$ with global minimum located at $`x=2.90`$; (B) CLM with $`\dot{x}=U(x)(xz)\lambda `$, $`\dot{z}=U(z)+(xz)+\lambda `$, $`\dot{\lambda }=xz`$, where $`x,z`$ (red and blue, respectively) are the two particle states with costate $`\lambda `$ (green). The CLM corresponds to the Lagrange programming network with cost function $`U(x)+U(z)`$ subject to $`x=z`$. Except when both initial states $`x(0),z(0)`$ are positive, the global minimum is always reached by this CLM. Similar results are obtained for other double well problems with broad and narrow minima.
Fig. 2 Supervised learning of multilayer perceptrons on a benchmark problem of training a sinusoidal function (green) given 20 training data points $`\{u_k,d_k\}`$ (blue circles) corrupted with Gaussian noise (zero mean, standard dev. 0.4) using a multilayer perceptron with one hidden layer consisting of 10 hidden units. Given the fact that the MLP consists of 30 unknown parameters, overfitting will occur with standard training methods, unless regularization or early stopping is applied. In (A) results from standard neural network training methods are shown: scaled conjugate gradient without early stopping leading to overfitting (red); best result by Bayesian learning with regularization of the cost function and automatic hyperparameter selection, leading to good generalization (blue). (B) shows that a CLM which optimizes a sum squared error on the training data without regularization of the cost function can produce a result (blue) which is comparable to the quality of the Bayesian learning solution in (A), provided the $`\sigma `$ value in $`p(𝐱(0))\mathrm{exp}(\frac{1}{2\sigma ^2}𝐱(0)^T𝐱(0))`$ is well chosen (here $`\sigma =0.1`$). (C) shows a result from a CLM for a bad choice of $`\sigma =5`$. The large variance among the resulting local minimizers of the CLM is clearly visible, in comparison with (B) where the variance is small with good generalization. In (D) the CLM evolution with $`q=20`$ of the sum squared error cost function is shown during optimization, related to the results of (B).
Fig. 3 CLM application to potential energy surface minimization of Lennard-Jones clusters: (A) CLM evolution of $`U_{LJ}^\delta `$ for $`q=50`$ on (LJ)38; (B) Importance of the choice of the initial distribution: $`\sigma =1`$ (red), $`\sigma =0.1`$ (blue) with post-processing by means of a quasi-Newton method resulting into $`\sigma =1`$ (magenta), $`\sigma =0.1`$ (green), respectively. Shown are the mean (dashed line), min and max (solid lines) over 7 CLM runs after sorting the energy values. The variety of solutions includes the global minimum -173.9284, visualized in (C). In (D) a CLM result with $`q=100`$ is shown for a larger cluster (LJ)150.

A

B

> Fig. 1

A

B

C

D

> Fig. 2

A

B

C

D

> Fig. 3


text:
# Orbital Dimerization and Dynamic Jahn-Teller Effect in NaTiSi2O6

## Abstract

We study with Raman scattering technique two types of phase transitions in the pyroxene compounds NaMSi2O6 (with M=Ti, V, and Cr). In the quasi one-dimensional S=1/2 system NaTiSi2O6 we observe anomalous high-temperature phonon broadening and large changes of the phonon energies and line-widths across the phase transition at 210 K. The phonon anomalies originate from an orbital order-disorder phase transition and these results –combined with theoretical considerations– indicate that the high temperature dynamical Jahn-Teller phase of NaTiSi2O6 exhibits a spontaneous breaking of translational symmetry into a dimerized, Jahn-Teller distorted, orbital ordered state under the formation of spin valence bonds. In S=1 NaVSi2O6 orbital degrees of freedom are strongly suppressed and the magnetic excitations are well described within a Heisenberg model, indicating that at TN=19K this system orders antiferromagnetically.

Introduction. Electrons in strongly correlated transition-metal compounds can be regarded as having separate spin, charge and orbital degrees of freedom. It is the interplay between these, combined with their coupling to the lattice, that gives rise to a wealth of possible spin, charge and orbital orderings, as observed for instance in many Colossal Magneto-Resistance manganites a1 , cuprates Kugel , titanates (e.g. LaTiO3 Keimer ) and vanadates (e.g. V2O3 Castellani  and LiVO2 Pen ).

Systems with orbital degeneracy are particularly interesting because orbitals couple to the lattice via the cooperative Jahn-Teller (JT) effect on one hand, and via superexchange interactions to the electronic spin on the other hand Kugel . Therefore at an orbital ordering phase transition the magnetic susceptibility and phonon properties will be affected at the same time. Experimentally, however, such an interrelation is seldom found: in general the dominant JT orbital-lattice coupling obscures the more subtle effects due to the superexchange.

Pyroxenes are large family of compounds (AMB2O6; A=alkali-metal, M=transition-metal, and B=Si, Ge) which structure consists of isolated quasi one-dimensional chains of edge-sharing MO6 octahedra, linked together by the corner-sharing BO4 tetrahedra, Fig. 1. Particularly, the compounds with active $`t_{2g}`$ orbitals (M=Ti,V, and Cr) are expected to have electronic interactions governed by both orbital degeneracy and anisotropy. Interestingly, the sodium-silicon system with titanium, NaTiSi2O6, is rather different from other pyroxenes as it lacks low-temperature AF order, and shows signs of the opening of a spin gap instead a2 ; a3 ; a30 ; a31 . Ti3+ corresponds to S=1/2, and since all Ti ions are in equivalent crystallographic positions, this material is a prime candidate to show a spin-Peierls (SP) phase transition. Indeed, its magnetic susceptibility a2  sharply decreases below 210 K, indicating a transition to a spin-singlet state, which could in principle be of the SP type. But already the high-temperature magnetic susceptibility data does not support such a SP scenario, since the phase transition occurs at a temperature that is higher then the maximum point of the Bonner-Fisher curve, which implies that the transition cannot solely be driven by magnetic fluctuations a2 . Based on structural analysis, Isobe et.al. anticipated that orbital dimerization might be responsible for such behavior a2 .

In this Letter, we report a Raman scattering study on a family of pyroxene compounds from which we conclude that orbital degrees of freedom play decisive role in physical properties of NaTiSi2O6, and that NaTiSi2O6 undergoes an orbital order phase transition at TOO=210 K. At high temperatures the fluctuations of the orbital degrees of freedom lead to a dynamic Jahn-Teller phase with anomalous phonon broadening and remnant antiferromagnetic (AF) spin fluctuations. We find that the dramatic drop of the magnetic susceptibility below TOO is accompanied with a structural change, just as one would expect for a canonical orbital ordering transition. These results, combined with the microscopic orbital-spin model that we derive, suggest that the quasi one-dimensional dynamical Jahn-Teller phase of NaTiSi2O6 exhibits a spontaneous breaking of the translational symmetry into an dimerized orbital ordered state with a spin gap due to the formation of spin valence bonds. The paper is organized as follows: first we present the Raman spectra and underline important observations, then we discuss the electronic structure of pyroxenes and derive a microscopic model. Finally, we demonstrate that all observed effects (including susceptibility data) can be accounted by our model.

Experiment.High-quality powder samples of NaMSi2O6 (M=Ti, V, and Cr) were prepared by a solid-state reaction of mixtures with an appropriate molar ratio of Na2MSi4O11, M and MO2. Details of the sample preparation are published elsewhere a2 . Raman spectra are measured in the backscattering configuration, using 514.5 nm line of an Ar-ion laser, micro-Raman system with a DILOR triple monochromator, and a liquid nitrogen cooled charge-coupled device detector. For the low-temperature measurements we used the Oxford continuous-flow cryostat with a 0.5 mm thick window. The focusing of the laser beam is realized with a long distance (10 mm focal length) microscope objective (magnification 50$`\times `$).

The pyroxenes crystallize in a monoclinic unit cell with the space group $`C2/c`$ a4 . The unit cell consists of four formula units (Z=4) with 40 atoms in all. The site symmetry of Na, M, Si, O1, O2 and O3 atoms are (4e), (4e), (8f), (8f), (8f) and (8f), respectively. Thus, the factor-group-analysis (FGA) yields: (Na, M)($`C_2`$) $`\mathrm{\Gamma }`$ = $`A_g`$ \+ $`A_u`$ \+ 2$`B_{1g}`$ \+ 2$`B_u`$; (Si, O1, O2, O3)($`C_1`$) $`\mathrm{\Gamma }`$ = 3$`A_g`$ \+ 3$`A_u`$ \+ 3$`B_g`$ \+ 3$`B_u`$. Summarizing these representations and subtracting the acoustic modes ($`A_u`$ \+ 2$`B_u`$), we obtain the following irreducible representations of NaMSi2O6 vibrational modes:

$`\mathrm{\Gamma }_{NaMSi_2O_6}^{opt.}`$ $`=`$ $`14A_g(xx,yy,zz,xz)+16B_g(xy,yz)`$
$`+`$ $`13A_u(𝐄||𝐲)+14B_u(𝐄||𝐱,𝐄||𝐳)`$
The unpolarized Raman spectra of NaMSi2O6 (M=Ti, V, Cr) are shown in Fig. 2. At the room temperature we find around 30 phonon modes as predicted by FGA. The Raman spectra of different compounds are similar as expected for isostructural materials, and the phonon modes can be crudely grouped into two energy regions. The modes in the spectral range below 500 cm-1 originate from the bond bending vibrations, whereas the higher frequency modes originate from the stretching vibrations. The highest energy modes are mainly due to the non-bridging Si-O ion vibrations, because of the shortest Si-O tetrahedral bonds. The modes of NaMSi2O6 at 1055/1032 cm-1 (Cr), 1042/1025 cm-1 (V), and 1042/1025 cm-1 (Ti), represent Si-O2 (see Fig. 2) antisymmetric/symmetric bond stretching vibrations, respectively. Their frequency difference scales as R-3 in a full accordance with difference between Si-O2 bond lengths in these materials. Similarly, the modes at 990/967 cm-1 (Cr), 972/954 cm-1 (V), and 965 cm-1(Ti), we assign as antisymmetric/symmetric pairs of Si-O1 bond stretching modes. This simplified mode assignment becomes inapplicable at lower frequencies due to more complicated normal coordinates of corresponding vibrations. However, besides similarities, we also observe two very important effects. First, most of the phonon line widths are dramatically increasing in the components with smaller spin values ($`V^{3+}S=1,Cr^{3+}S=3/2`$). Second, the spectra of NaTiSi2O6, due to the large phonon broadening, show effectively less phonon modes then expected by FGA, and observed in other pyroxenes. The latter effect may also be regarded as a consequence of a high temperature ”higher symmetry lattice state” of NaTiSi2O6.

By lowering the temperature we find a dramatic change in the phonon Raman spectra of NaTiSi2O6 around 210 K, see Fig. 3a,b. The modes at 221 and 209 cm-1 exhibit the anti-crossing behavior, see inset Fig.3a, which indicates that these two modes belong to the same symmetry class below $`T_c`$. A similar effect is observed for all other such pair of modes (one more pair is also visible in Fig. 3a at about 180 and 190 cm-1) in the spectra, which suggests that in the low-T phase of NaTiSi2O6 all Raman active modes are of the same (Ag) symmetry type. Furthermore, our reflectance measurements (will be published separately, Ref. a6 ) show that the infrared active phonon modes can be distinguished from the Raman modes, since the center of inversion remains to be the symmetry element in the low-T phase. Accordingly, we uniquely determine the space group of the low-T phase of NaTiSi2O6 to be P$`\overline{1}`$. This space group requires the translation symmetry breaking of the TiO6 chain (the tentative distortion pattern of the low-T phase is shown in Fig. 1) which is in agreement with preliminary X-ray diffraction and neutron scattering data a7 .

Moreover, we show in details the temperature dependence of the structure around 970 cm-1, Fig. 3b. The mode at about 946 cm-1 softens by about 10 cm-1, while mode at 966 cm-1 ”splits”, and hardens by 25 cm-1, see inset in Fig. 3b. The full width at half maximum (FWHM) of the 946 cm-1 (circles) and 966 cm-1 (squares) phonons (the FWHM is estimated from Lorenzian-fit) increase up to the maximum value at about 210 K, and then decreases to the saturation value which is much smaller then the T=300 K value, see right inset of Fig. 3b. In fact, all Raman-active phonons exhibit anomalies at the temperature which coincides with $`T_c`$ obtained from susceptibility measurements a2  indicating that the magnetic ordering is accompanied with the structural phase transition just as one would expect for an canonical orbital-ordering transition. Besides, behavior of FWHM implies that the bond fluctuations are considerably larger in the high-T then in low-T phase(proximity of the structural phase transition induces the largest fluctuations, producing the maximum FWHM at TOO). Due to the JT effect bond fluctuations emphasize the strong orbital character of disorder above the phase transition temperaturea771 ; a772 .

Contrary to S=1/2 compound, in the low temperature Raman spectra of NaVSi2O6, see Fig. 3c, we did not observe the phonon anomalies. Instead, at low temperatures, broad asymmetric features, typical for the two-magnon excitations in Heisenberg AF ($`H=_{i,j}J_{i,j}𝐒_i𝐒_j,J_{i,j}>0`$), appear around 260 and 400 cm-1. According to the linear Anderson’s approximation, differential cross section for the two-magnon (tm) Raman scattering is proportional to the two-magnon density of states (DOS), $`I_{tm}_𝐤\delta (\omega 2\omega (𝐤))`$. The dispersion relation is obtained assuming two-dimensional magnetic structure with different exchange interactions along the chains (J), and perpendicular (J) to the chains: $`\omega (𝐤)=2S\sqrt{(J_{}+J_{})^2(J_{}cos(ka)+J_{}cos(kb))^2}`$. The calculated two-magnon DOS is shown in top panel of Fig. 3c. Two singularities at energies of about 260 and 390 cm-1 are obtained for J$`{}_{}{}^{}=85`$ cm-1 and J$`{}_{}{}^{}=13`$ cm-1 (120 K and 18 K, respectively), in a very good agreement with experiment. The maximum of the susceptibility curve in NaVSi2O6 is at about 100 K a7  (giving $`J_{}`$ 80 K a301 ) which is in good agreement with our estimation of the J. Thus, we conclude that NaVSi2O6 exhibits an AF ordering. Similar AF phase has been already observed in other S=1 chain compounds LiVGe2O6 a3 , and NaVGe2O6 a31 .

Pyroxenes: theoretical. Next we discuss the electronic structure. For a single octahedron the cubic crystal field splits the Ti,V or Cr 3d-states into low lying $`t_{2g}`$ states, with 1,2 or 3 electrons, respectively, and empty states of $`e_g`$ symmetry at higher energy. The low energy electronic properties are governed by the three-fold degenerate $`t_{2g}`$ states (the relevant states are $`|xy`$, $`|yz`$ and $`|zx`$a701 .

The Coulomb interaction $`U`$ between electrons on the same transition metal atoms is large, so the exchange interactions can be determined by a second order perturbation expansion in the electron hopping parameters. The problem is further reduced by considering the symmetry allowed hopping paths in the chain geometry, schematically represented in Fig. 4. If we consider orbitals on two sites in the same $`xy`$ plane, then only the hopping between $`|xy`$ orbitals is relevant (see Fig. 4a). For sites in the $`yz`$ plane, the $`|yz`$ orbitals are relevant, Fig. 4b. In the present geometry (note that x, y, and z in Fig. 4 do not correspond to a, b, and c axis in Fig. 1) there are no transition metal atoms in the same chain that are also in the same $`xz`$ plane. The $`|xz`$ orbitals are therefore non-bonding and can be considered inert on this level of approximation (a tightbinding parametrization a5a  shows that other overlap integrals either vanish by symmetry or are more than factor five smaller), see Fig. 4c. For the S=1/2 Titanium system we then obtain the Hamiltonian a5b ,

$$H^{\mathrm{Ti}}=\mathrm{J}^{\mathrm{Ti}}\underset{i,j}{}𝐒_i𝐒_j[\frac{1}{4}+T_i^zT_j^z+\frac{(1)^i}{2}(T_i^z+T_j^z)],$$
(1)
where we use the orbital operators $`T`$ ($`T_i^z=1/2`$ corresponds to an occupied $`|xy`$ orbital and $`T_i^z=1/2`$ to an occupied $`|yz`$ orbital on site $`i`$), and $`i,j`$ are neighboring sites. For the S=1 Vanadium and S=3/2 Chromium system the orbital degree of freedom vanishes ($`T=0`$) as the both of the active orbitals are occupied, so for the spin system we obtain a simple Heisenberg Hamiltonian $`H^{\mathrm{V}/\mathrm{Cr}}=\mathrm{J}^{\mathrm{V}/\mathrm{Cr}}_{i,j}𝐒_i𝐒_j`$. This immediately explains the uniqueness of the Ti (S=1/2) system, and the antiferromagnetism in NaVSi2O6 (in fact in all S $`1/2`$ pyroxenes).

Interpretation of experiments. The ground state of Hamiltonian (1), is clearly a ferro-orbital state, with spin-singlets on each bond, where the energy per dimer is $`3\mathrm{J}^{\mathrm{Ti}}/4`$. The state with all $`|xy`$ occupied is degenerate with the state with all $`|yz`$ occupied, see Fig. 4. Those states do differ, however, because the dimerization pattern along the chain is shifted by one lattice spacing. At zero temperature the system is condensed in either one of these two dimerized orbital ordered states and the translation symmetry is broken. This explains the structural change at $`T_c`$ a7 , the symmetry change and the energy shifts of the phonon excitations in the Raman spectra of NaTiSi2O6, and the observation of a large spin gap in the susceptibility measurements a2 .

At high temperatures, due to strong orbital-lattice coupling and JT effect, the orbital fluctuations produce dramatic effects. Above TOO the JT distortion disappears: the high-T orbital disordered phase of NaTiSi2O6 may be regarded as an orbital fluctuating phase, and on larger time scales the TiO6 octahedra appear to be undistorted. This causes the crystal to be effectively more symmetric, which is in agreement with the observation of less-then-expected phonon modes in the room temperature Raman spectra. As a consequence, orbital fluctuations induce large phonon broadenings (of course, the modes with Ti-O bonds in their normal coordinates will be mostly affected), as they are indeed found in the room temperature Raman spectra of NaTiSi2O6, see Fig. 2. In that respect the high-T phase of NaTiSi2O6 resembles a dynamical JT phase (where the phonon broadening is a signature of the melted static lattice distortions a771 ; a772 ). The spin fluctuations above the orbital ordering temperature a2  can be easily understood in a straightforward mean-field approximation of Hamiltonian (1), where we can decouple orbitals and spins. If the long range order is absent in the orbital sector, we have for the expectation values $`T_i^z=0`$, and $`|T_i^zT_j^z|<1/4`$, so that the effective exchange constants –and therefore the fluctuations in the spin-sector– are still antiferromagnetic. Thus, the nature of the phase transition in NaTiSi2O6 can be established: it corresponds to an orbital order-disorder phase transition with appropriate concomitant magnetic and lattice changes.

In the case of two electrons per site (NaVSi2O6), the orbital fluctuations are strongly suppressed due to inert property of $`|xz`$ orbitals, and the corresponding phonon broadenings are much smaller, see Fig. 2. NaCrSi2O6 has a fully polarized t2g core, no orbital degrees of freedom, and no anomalous phonon broadenings in the Raman spectra.

Conclusions. We report a study of the two types of the phase transitions, observed in pyroxene family, by analyzing the Raman-active phonon and magnon excitations, and their temperature dependence. We find that spin S=1/2 compound, NaTiSi2O6, exhibits a transition at about 210 K that we assign to be an orbital order-disorder phase transition. It originates from the instability of the high temperature orbital fluctuating - dynamical Jahn-Teller phase, towards a dimerized orbital ordered state, which is accompanied by a lattice distortion and by spin valence bond formation. The spin S=1 system, NaVSi2O6, on the contrary, does not show this type of instability and orders as a Nèel-type antiferromagnet below TN = 19 K in agreement with our microscopic spin-orbital model for the pyroxenes.

M.J.K. and Z.V.P. acknowledge support from the Research Council of the K.U. Leuven and DWTC. The work at the K.U. Leuven is supported by Belgian IUAP and Flemish FWO and GOA Programs.

\*milan.konstantinovic@fys.kuleuven.ac.be

\**Present address: Materials Science Institut, University of Valencia, Poligono La Coma, 46980 Paterna (Valencia), Spain


text:
# The Role of Mass and External Field on the Fermionic Casimir Effect

## I Introduction

As well as the Lamb-Retherford shift of atomic energy levels and the electron’s anomalous magnetic moment, the Casimir force between two parallel perfectly conducting plates is one of the most remarkable manifastations of quantum vacuum fluctuantions. It was first predicted on theoretical grounds by H. B. G. Casimir in 1948  and experimentally verified on a qualitative level by Sparnaay  ten years later. Recentely, high accuracy experiments have been performed by Lamoreaux  and by Mohideen and Roy . For a more detailed account on the subject there are excellent reviews in the literature .

In 1975, employing Casimir essential ideas, Johnson  investigated the effects of boundaries on a massless Dirac field in the context of MIT-bag model  and found an energy density shift of the same order of magnitude as that obtained by Casimir for the electromagnetic field. Addopting Johnson’s extend approach to Casimir effect, in order to allow for other quantum fields, many authors have investigated the effects of different boundary conditions on the fields considered . Hence, one can say that a modern view of Casimir effect might take into account those effects caused by non-trivial space topologies on the vacuum of quantum fields.

It is worthnoting that in this general context some Casimir setups (field + boundary condition \+ external sources) present quite complicated final expressions for the Casimir energy density, which ultimately obscure any possible physical interpretation. In particular, the results obtained in  for the massive spinor field indicate a mass dependent energy density which calls for a deep investigation. Similar difficulties also appear in the case considered in  where the fermion field is also subjected to an external magnetic field.

In order to handle the above mentioned shortcomings, we present here an alternative treatment which allows us to extract new information concerning the role of mass and that of an external field to the fermionic Casimir effect. This is achived by a suitable combination of the method of analytic regularization in the context of gamma function representation (also called $`\alpha `$-representation ) and the use of the well-known Euler-Mclaurin summation formula .

The main ideas of our construct are introduced in the next section, where we consider the well-stablished electromagnetic Casimir effect. In section III, the same procedure is employed to the case of a massive spinor field with boundaries analogous to that considered in  and . In section IV, the effects of an external constant and homogeneous magnetic field is also considered and the connection with the so-called Euler-Kockel-Heisenberg Effective Lagrangian density  is addressed. Finaly, in section V, we make some concluding remarks pointing out directions of future investigations.

## II The Case of Electromagnetic Fields

As was originally considered by Casimir in 1948, the divergent vacuum energy density associated with the electromagnetic field is given by

$$\epsilon _0=\frac{1}{8\pi ^2L}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}𝑑k_x𝑑k_y\frac{1}{[k_x^2+k_y^2+(an)^2]^{1/2}},a=\pi /L$$
(1)
where the discrete values of $`k_z`$ reflects the boundary condition imposed by two parallel conducting plates with area $`A=l^2`$ and separated by a distance $`L`$ in such a way that $`l>>L`$ ($`a=\pi /L`$) .

In (1), the integrals are quadraticaly divegent quantities which claim for a consistent regulariazation prescription. Despite the familiar regularization methods found in the literature  , we shall consider here a quite different one. It consists in the conbination of the analytic regulariazation scheme, using the gamma function integral representation, with the so-called Euler-Mclaurin summation formula . To see how this works, we start by taking the analytic extension of the integrand in (1), which turns out to be a regular functional. This is achieved by means of the gamma function integral representation

$$\frac{1}{A^{1+\delta }}=\frac{1}{\mathrm{\Gamma }(1+\delta )}_{0^+}^{\mathrm{}}𝑑\eta \eta ^\delta \mathrm{e}^{A\eta },$$
(2)
valid for $`\delta >1`$, which allows us to rewrite (1) as

$$(\epsilon _0)^\mathrm{R}=\frac{1}{2(2\pi )^2L}\frac{\pi }{\mathrm{\Gamma }(1/2+\delta )}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}_{0^+}^{\mathrm{}}𝑑\eta \eta ^{5/2+\delta }\mathrm{e}^{(an)^2\eta },$$
(3)
where the gaussian integrals in $`k_x`$ and $`k_y`$ have already been calculated. Note that, for $`\delta 0`$, we expect to recover the original theory, i.e., expression (1). This will be done only at the end of the calculations.

The divergent sum over $`n`$ appearing in (3) is performed by means of the Euler-Mclaurin summation formula

$`{\displaystyle \underset{n=M}{\overset{N}{}}}f(n)`$ $`=`$ $`{\displaystyle _M^N}F(x)𝑑x+{\displaystyle \frac{1}{2}}[f(N)+f(M)]+{\displaystyle \underset{k=1}{\overset{K}{}}}{\displaystyle \frac{B_{2k}}{2k!}}[F^{2k1}(N)F^{2k1}(M)]`$ (5)
$`+{\displaystyle \frac{1}{(2K+1)!}}{\displaystyle _N^M}B_{2K+1}(x[x])F^{2K+1}(x)𝑑x`$
where $`B_mB_m(0)`$ and the $`B_m(x)`$ are the Bernoulli polynomyals. We preserve here the same notation used in . The last term in (4), also called the remainder term, vanishes if $`F(z)`$ is an entire function. Further, if $`n`$ is integer and $`MnN`$, then $`F(n)=f(n)`$. In the present context, identifying the entire function $`f(n)`$ with

$$f(n)=\mathrm{e}^{(an)^2\eta }$$
(6)
and, since $`0n\mathrm{}`$, we are allowed to rewrite (3) as

$`(\epsilon _0)^\mathrm{R}`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2L}}{\displaystyle \frac{\pi }{\mathrm{\Gamma }(1/2+\delta )}}\{{\displaystyle _{0^+}^{\mathrm{}}}d\eta \eta ^{5/2+\delta }[{\displaystyle _0^{\mathrm{}}}dnf(n){\displaystyle \frac{1}{2}}f(n)|_n\mathrm{}`$ (8)
$`+{\displaystyle \frac{1}{12}}f^{}(n)|_n\mathrm{}{\displaystyle \frac{1}{12}}f^{}(n)|_{n0}{\displaystyle \frac{1}{720}}f^{\prime \prime \prime }(n)|_n\mathrm{}+{\displaystyle \frac{1}{720}}f^{\prime \prime \prime }(n)|_{n0}+\mathrm{}]\}`$
where $`f^{}(n)`$ means the first derivative of (5) with respect to $`n`$ and so on. Calculating the derivatives and taking the corresponding limits we see that the only non-vanishing terms give

$`(\epsilon _0)^\mathrm{R}`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2L}}{\displaystyle \frac{\pi }{\mathrm{\Gamma }(1/2+\delta )}}\{{\displaystyle _{0^+}^{\mathrm{}}}d\eta \eta ^{3+\delta }{\displaystyle \frac{\sqrt{\pi }}{2a}}+{\displaystyle \frac{1}{12}}\left({\displaystyle \frac{\mathrm{\Gamma }(1/2+\delta )}{(a^2n^2)^{1/2+\delta }}}(2a^2n)\right)_n\mathrm{}`$ (10)
$`{\displaystyle \frac{1}{720}}\left({\displaystyle \frac{\mathrm{\Gamma }(3/2+\delta )}{(a^2n^2)^{3/2+\delta }}}(2a^2n)^3\right)_n\mathrm{}{\displaystyle \frac{3}{720}}\left({\displaystyle \frac{\mathrm{\Gamma }(1/2+\delta )}{(a^2n^2)^{1/2+\delta }}}(2a^2)^2n\right)_n\mathrm{}\}`$
where (2) was used. In going from expression (6) to (7) we notice that only the ($`n\mathrm{}`$)- terms contribute to the energy density. Usually, the methods found in the literature extract information arising from the ($`n0`$)-terms and, as we will show in the next section, this generates quite different final results. However, we must stress that, while the exponential function in (5) is an analytical function over the entire complex plane, the power function in the integrand of (1) is a multiple-valued function, which has a branch cut along the real axis .

Now, in order to obtain a consistent result with the original theory we now take the limit $`\delta 0`$ in (7). First, we must handle the divergent contribution arising from the second term in the curly brackets. Such divergence is eliminated using the freedom in the choice of $`\delta `$ in (7). In fact, for consistency with (3), $`\delta `$ are constrained to be greater than $`1/2`$. However, to obtain the correct final result the considered region in the complex plane must be analytic continued to $`\delta 1`$. Hence, as the limit $`n\mathrm{}`$ is performed, the second term contribution turns out to be zero. On the other hand, the contribution from the remaining terms, which becomes $`n`$-independent when the limit $`\delta 0`$ is taken, gives

$$\epsilon _0=\frac{1}{4(2\pi )^2}_0^{\mathrm{}}𝑑\eta \eta ^3+\frac{\pi ^2}{720L^4}.$$
(11)
Finally, calculating the corresponding energy for the whole space (which is equivalent to taking $`L\mathrm{}`$ in the above expression) and subtracting from it expression (8), yields

$`\mathrm{\Delta }\epsilon `$ $`=`$ $`(\epsilon _0)_L\mathrm{}\epsilon _0`$ (12)
$`=`$ $`{\displaystyle \frac{\pi ^2}{720L^4}}`$ (13)
which is the expected vacuum energy density derived by Casimir in 1948. Since there is no a priori reason for assuming that the vacuum energy in the presence of boundaries is greater than in their absense, we have defined $`\mathrm{\Delta }\epsilon `$ in such a manner that (9) indicates that the force between the plates is in fact attractive ***Once we adopt this definition, it must be preserved if we consider other topologies and boundary conditions, in order to be able to distinguish attractive from repulsive effects., as has been confirmed by precise experiments  \- .

## III The Massive Spinor field: $`S^1\times R^3`$

The case of noninteracting spinor fields subjected to $`S^1\times R^3`$ space topology with twisted and untwisted spin connections was first considered by DeWitt, Hart and Ishan in 1979  and one year later it was generalized to the massive case by Ford , who obtained an intrincate mass dependent expression for the vacuun energy density. In this section we intend to rederive these results using the procedure developed in the last section. As will be shown, the effect of the fermion mass to the Casimir energy density is, in fact, null.

### A Untwisted Case

For an untwisted spinor field, the vacuun energy density is given by

$$\epsilon _0^{\mathrm{unt}}=\frac{1}{(2\pi )^2L}\underset{\mathrm{}}{\overset{\mathrm{}}{}}_{\mathrm{}}^{\mathrm{}}𝑑p_x𝑑p_y[m^2+p_x^2+p_y^2+p_z^2]^{1/2}.$$
(14)
where $`p_z=a^2n^2`$, with $`n=0,\pm 1,\pm 2,\mathrm{}`$ and $`a=2\pi /L`$. As in the electromagnetic case, the above quantity is divergent and must also be regularized. Using the same procedure employed in the last section, we write (10) as

$`(\epsilon _0^{\mathrm{unt}})^R`$ $`=`$ $`{\displaystyle \frac{1}{2\pi L}}{\displaystyle \frac{1}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{5/2+\delta }\mathrm{e}^{(m^2+a^2n^2)\eta }`$ (15)
$`=`$ $`{\displaystyle \frac{1}{2\pi L}}{\displaystyle \frac{1}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{5/2+\delta }\left[\mathrm{e}^{m^2\eta }+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{e}^{(m^2+a^2n^2)\eta }\right],`$ (16)
where we have already performed two gaussian integrals in $`p_x`$ and $`p_y`$. Instead of use the Abel-Plana formula - we may invoke the Euler-Mclaurin summation formula (4), in order to find an expression similar to (6), where the only difference rely upon the definition of the entire function $`f(n)`$. Here

$$f(n)=\mathrm{e}^{(m^2+a^2n^2)\eta }.$$
(17)
Following the same steps done in the last section, it is straightforward to obtain

$`(\epsilon _0^{\mathrm{unt}})^R`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )L}}{\displaystyle \frac{2}{\mathrm{\Gamma }(1/2+\delta )}}\{{\displaystyle _{0^+}^{\mathrm{}}}d\eta \eta ^{3+\delta }\mathrm{e}^{m^2\eta }{\displaystyle \frac{\sqrt{\pi }}{2a}}+{\displaystyle \frac{1}{12}}\left({\displaystyle \frac{\mathrm{\Gamma }(1/2+\delta )}{(m^2+a^2n^2)^{1/2+\delta }}}(2a^2n)\right)_n\mathrm{}`$ (18)
$``$ $`{\displaystyle \frac{1}{720}}\left({\displaystyle \frac{\mathrm{\Gamma }(3/2+\delta )}{(m^2+a^2n^2)^{3/2+\delta }}}(2a^2n)^3\right)_n\mathrm{}{\displaystyle \frac{3}{720}}\left({\displaystyle \frac{\mathrm{\Gamma }(1/2+\delta )}{(m^2+a^2n^2)^{1/2+\delta }}}(2a^2)^2n\right)_n\mathrm{}\}.`$ (19)
Analyzing the structure of (13), we see that the fermion mass appearing in the denominator of the last three terms may be neglected since in those terms $`n\mathrm{}`$. This point is crucial since the limit accounts for the partial elimination of the fermion mass. In fact, $`m`$ remains in the kernel of the first term but, as will shown later, this term will be cancelled against the Minkowiski vacuum energy . The divergence of the second term in (7) also appears and we should say that it is a common feature of our approach. Again, it is eliminated by taking an analytic extension analogous to that employed in the electromagnetic case. As a final result we have

$$\epsilon _0^{\mathrm{unt}}=\frac{1}{8\pi ^2}_0^{\mathrm{}}𝑑\eta \eta ^3\mathrm{e}^{m^2\eta }\frac{\pi ^2}{720L^4}.$$
(20)
Subtracting (14) from the corresponding usual Minkowiski vacuun energy (which correspond to taking $`L\mathrm{}`$ in (14)) we find the fermionic Casimir energy density

$$\mathrm{\Delta }\epsilon ^{\mathrm{unt}}=\frac{2\pi ^2}{45L^4},$$
(21)
which is in complete agreement with , the only difference being a factor $`2`$, which reflects the four-component spinor field representation we are using. De Witt et al considered a two-component spinor field.

### B Twisted Case

Expression (14) clearly shows the independence of the Casimir energy density with respect to the fermion mass. This feature also occurs when twisted boundary condition is considered. In this case $`p_z=(2n+1)a`$, with $`n=0,\pm 1,\pm 2,\mathrm{}`$ and $`a=\pi /L`$. The regulated vacuun energy density is now given by

$$(\epsilon _0^{\mathrm{twi}})^R=\frac{\sqrt{\pi }}{(2\pi )^2L}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}_0^{\mathrm{}}𝑑\eta \eta ^{5/2}\mathrm{e}^{(m^2+(2n+1)^2a^2)\eta }.$$
(22)
Since the above integral is regular, we are allowed to interchange the sum and the integral, and than use the following mathematical trick

$$\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{e}^{(m^2+(2n+1)^2a^2)\eta }=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{e}^{(m^2+n^2a^2)\eta }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{e}^{(m^2+(2n)^2a^2)\eta }.$$
(23)
In this way the problem of solving the twisted case reduces to that of computing two terms proportional to that in the untwisted case. In fact, we have

$`\epsilon _0^{\mathrm{twi}}`$ $`=`$ $`{\displaystyle \frac{1}{2^3}}\epsilon _0^{\mathrm{unt}}\epsilon _0^{\mathrm{unt}}`$ (24)
$`=`$ $`{\displaystyle \frac{7}{8(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^3+2{\displaystyle \frac{7\pi ^2}{360L^4}}.`$ (25)
Again, subtracting this result from that where $`L\mathrm{}`$, we obtain

$$\mathrm{\Delta }\epsilon _0^{\mathrm{twi}}=2\frac{7\pi ^2}{360L^4},$$
(26)
which coincide with the familiar result found in the literature for the massless fermionic Casimir effect .

## IV The Influence of an External Field

Another interesting problem to be analysed using the present construct is related to influence of an external magnetic field on the Casimir energy associated to the Dirac field. This problem has been recentely proposed in the context of Effective Quantum Electrodynamics using the so-called Schwinger proper-time method . However, a clear answser concerning the role of the external field on the fermionic Casimir energy density is yet an open problem which deserves further investigation. The purpose of this section is to implement, in the context of Weisskopf method  , the prescription presented in the previous sections in order to get a better understanding of the above mentioned problem.

We restrict our calculation to the case where an untwisted massless spinor field is subjected to an external constant uniform magnetic field. As is well known , the negative energy levels for an electron of charge $`e=|e|`$ in the presence of an uniform and constant magnetic field $`\mathrm{H}_z=\mathrm{H}`$ is given by

$$ϵ_{𝐩,\sigma }^{()}=\sqrt{m^2+(2n+1\sigma )|e|\mathrm{H}+p_z^2},$$
(27)
where $`n=0,1,2,3\mathrm{}`$ and $`\sigma =\pm 1`$. Taking into account the density of states in the interval $`dp_z`$

$$\frac{|e|\mathrm{H}}{2\pi }\frac{dp_z}{2\pi }$$
(28)
and the fact that all the levels except $`n=0,\sigma =1`$ are doubly degenerate (the levels $`n,\sigma =1`$ and $`n+1,\sigma =1`$ coincide), we obtain the energy density of vacuum electrons,

$`\epsilon _0^{}`$ $`=`$ $`{\displaystyle \underset{𝐩,\sigma }{}}ϵ_{𝐩,\sigma }^{()}`$ (29)
$`=`$ $`{\displaystyle \frac{|e|\mathrm{H}}{(2\pi )^2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}\left\{\sqrt{p_z^2}+2{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\sqrt{2|e|\mathrm{H}n+p_z^2}\right\}𝑑p_z.`$ (30)
where $`p_z`$ turned out to be a discrete quantity in virtue of the untwisted $`S^1\times R^3`$ space topology we are assuming. Using (2), the energy density (22) may be rewritten in the more convenient form,

$`(\epsilon _0^{})^R`$ $`=`$ $`{\displaystyle \frac{|e|H}{2\pi L}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\eta \eta ^{3/2+\delta }\left({\displaystyle \underset{n^{}=\mathrm{}}{\overset{\mathrm{}}{}}}\mathrm{e}^{(2|e|\mathrm{H}n+a^2n^2)\eta }\right)`$ (31)
$`=`$ $`{\displaystyle \frac{|e|H}{2\pi L}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\eta \eta ^{3/2+\delta }\mathrm{e}^{\alpha \eta }\left(f(0)+2{\displaystyle \underset{n^{}=1}{\overset{\mathrm{}}{}}}f(n^{})\right)`$ (32)
with $`2|e|\mathrm{H}n=\alpha (n)\alpha `$,

$$f(n^{})=\mathrm{e}^{(an^{})^2)\eta },$$
(33)
and $`n^{}=0,\pm 1,\pm 2,\mathrm{}`$, $`a=2\pi /L`$. Aplying the Euler-Maclaurin formula (4) and performing the corresponding derivatives and limits, we arrive at

$`(\epsilon _0^{})^R`$ $`=`$ $`{\displaystyle \frac{\sqrt{\pi }}{(2\pi )^2}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{2+\delta }{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{e}^{\alpha \eta }`$ (37)
$`{\displaystyle \frac{2}{3L^3}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^{}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{1/2+\delta }\mathrm{e}^{(an)^2\eta }{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{e}^{\alpha \eta }\right\}_n^{}\mathrm{}`$
$`{\displaystyle \frac{2^4(2\pi )^5}{720L^7}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^3{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{3/2+\delta }\mathrm{e}^{(an)^2\eta }{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{e}^{\alpha \eta }\right\}_n^{}\mathrm{}`$
$`+{\displaystyle \frac{12(2\pi )^3}{720L^5}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^{}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{1/2+\delta }\mathrm{e}^{(an)^2\eta }{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\mathrm{e}^{\alpha \eta }\right\}_n^{}\mathrm{}.`$
The sum in the integrands can be eliminated by noting that

$$\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{e}^{\alpha \eta }=\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{e}^{2|e|\mathrm{H}\eta }=\mathrm{coth}(|e|\mathrm{H}\eta ).$$
(38)
Furthermore, assuming the weak field regime ($`\mathrm{H}<<1`$) we are allowed to expand the kernel of the integrals in (25), namely,

$`(\epsilon _0^{})^R`$ $`=`$ $`{\displaystyle \frac{\sqrt{\pi }}{(2\pi )^2}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{2+\delta }\mathrm{coth}(|e|\mathrm{H}\eta )`$ (42)
$`{\displaystyle \frac{2}{3L^3}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^{}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{1/2+\delta }\mathrm{e}^{(an)^2\eta }\left({\displaystyle \frac{1}{|e|\mathrm{H}\eta }}+{\displaystyle \frac{|e|\mathrm{H}\eta }{3}}+\mathrm{\Sigma }\right)\right\}_n^{}\mathrm{}`$
$`{\displaystyle \frac{2^4(2\pi )^5}{720L^7}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^3{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{3/2+\delta }\mathrm{e}^{(an)^2\eta }\left({\displaystyle \frac{1}{|e|\mathrm{H}\eta }}+{\displaystyle \frac{|e|\mathrm{H}\eta }{3}}+\mathrm{\Sigma }\right)\right\}_n^{}\mathrm{}`$
$`{\displaystyle \frac{12(2\pi )^3}{720L^5}}{\displaystyle \frac{|e|\mathrm{H}}{\mathrm{\Gamma }(1/2+\delta )}}\left\{n^{}{\displaystyle _0^{\mathrm{}}}𝑑\eta \eta ^{1/2+\delta }\mathrm{e}^{(an)^2\eta }\left({\displaystyle \frac{1}{|e|\mathrm{H}\eta }}+{\displaystyle \frac{|e|\mathrm{H}\eta }{3}}+\mathrm{\Sigma }\right)\right\}_n^{}\mathrm{}.`$
where

$$\mathrm{\Sigma }\underset{k=2}{\overset{\mathrm{}}{}}\frac{2^{2k}B_k}{(2k)!}(|e|\mathrm{H}\eta )^{2k1},$$
(43)
and the $`B_k`$’s are the Bernoulli numbers.

We are now in position to perform, term by term in the expansion of integral (27), the limit $`n^{}\mathrm{}`$. After a straightforward calculation we obtain

$$\epsilon _0^{}=\frac{|e|\mathrm{H}}{8\pi ^2}_0^{\mathrm{}}𝑑\eta \eta ^2\mathrm{coth}(|e|\mathrm{H}\eta )\frac{2\pi ^2}{45L^4},$$
(44)
where the same kind of analytic extension made in the previous sections was performed in the manipulation of the first term in the second line of (27). Again, it gives no contribution.

Finaly, the energy density of the empty space may be obtained by taking the limit of zero field and infinite volume in (29). We must subtract (29) from this quantity, obtaining

$$\mathrm{\Delta }\epsilon _0=\frac{1}{8\pi ^2}_{0^+}^{\mathrm{}}\frac{d\eta }{\eta ^3}\left\{|e|H\eta \mathrm{coth}(|e|H\eta )1\right\}+\frac{2\pi ^2}{45L^4},$$
(45)
which clearly shows the influence of the external magnetic field to the fermionic Casimir effect. It must be noted that the above expression recovers (15) in the limit of zero magnetic field. In addition, the first term in (30) might be recognized as the Euler-Kockel-Heisenberg correction to the effective Lagrangian density, which accounts for the nonlinear effects induced by the external field in effective quantum electrodynamics  \- . It provides exactly the same contribuction obtained when the limit $`L\mathrm{}`$ is considered, i.e., the contribution from the boundaries just add a field independent ammount to the E-K-H effective Lagrangian density. The independence of both effects clarify the physics governing the behaviour of quantum fields under the influence of external fields and/or boundaries conditions. The generalization of the above calculation to the twisted case is immediate as well as the inclusion of the fermion mass.

## V Concluding Remarks

Using an approach based on the combination of analytic regularization method throught $`\alpha `$-representation and the Euler-Maclaurin summation formula we had rederived the electromagnetic and the fermionic Casimir energy densities. The later, which comprises the main results of the present article, was considered in the case of $`S^1\times R^3`$ space topology, where the role played by the fermion mass and that of an external field on the Casimir energy density were fully investigated.

As was shown in section III, the present approach provided a powerfull way to dealwith in each step of the calculation, the divergences inherent to the theory. It was found that the fermion mass doesn’t play any influence on the twisted and untwisted fermionic Casimir energy densities, which is in contrast with the first results obtained by Ford . Experiment may provide the final answer.

We have also seen that, when an external magnetic field is considered, its effect on the Casimir energy density appears as an $`L`$-independent term (which was ultimately identifyed with the well known Euler-Kockel-Heisenberg Effective Lagrangian density) plus a term identical to that obtained when the external field is absent, expresion (15). This result clearly shows the independece of the external field on the boundary conditions, although this seems to be in apparent disagreement with the results found in .

Finally, it must be emphasized that the present construct is a simple and easily generalizable method to reexamine many other phenomena. Among these are those related to the Effective Quantum Electrodynamics in the context of the “old fashioned” Weisskopf’s method , recentely readdressed .

Acknowledgements. JLT acknowledges IFT/UNESP for the hospitality. LCC is grateful to FAPESP for the financial support.


text:
# “Weather” Variability Of Close-in Extrasolar Giant Planets

## 1. Introduction

The focus of extrasolar planet research has broadened to now include the characterization of their physical properties, as shown by the recent sodium detection in the atmosphere of HD 209458 b (Charbonneau et al. 2002). Atmospheric circulation is expected to play a key role in determining a number of observational characteristics of EGPs, including their albedo and transmission spectrum (see, e.g., Seager & Sasselov 1998, 2000; Sudarsky et al. 2000; Brown 2001). This is especially true for close-in EGPs, which are thought to be tidally-locked to their parent star and irradiated on one side only: circulation will be essential in redistributing heat from the day to the night side on these planets, thus determining to a large extent how they will appear to the distant observer (Cho et al. 2002a,b; Showman & Guillot 2002).

Recently, we have presented a set of detailed shallow-water numerical simulations of the atmospheric flow on HD 209458 b (Cho et al. 2002a,b), currently the only EGP with known mass and radius from the transit light curves and radial velocity measurements (Charbonneau et al. 2000; Henry et al. 2000; Mazeh et al. 2000; Jha et al. 2000; Brown et al. 2001). These simulations suggest that, contrary to the simple day/night (hot/cold) picture, the circulation on this planet is characterized by two moving circumpolar vortices and a small number of latitudinal bands/jets. The vortices act as dynamically distinct thermal spots whose motion around the poles generates variability as seen by an observer interested in quantities integrated over the planetary disk (or circumference).

It is possible to determine the general features of the circulation pattern expected within the framework of shallow-water dynamics by specifying the two dimensionless numbers – Rossby ($`R_o`$) and Burger ($`B_u`$) – for the atmospheric flow. In this Letter, we estimate a range of likely $`R_o`$ and $`B_u`$ values for known close-in EGPs and conclude that their atmospheric circulation pattern should be qualitatively similar to that of HD 209458 b. In §2, we recall how the atmospheric flow pattern can be characterized by the knowledge of $`R_o`$ and $`B_u`$. In §3, we describe the sample of close-in EGPs selected for our study and how we estimate likely values for various global planetary parameters entering into the definition of $`R_o`$ and $`B_u`$. Finally, our results and conclusions are presented in §4.

## 2. Turbulent Shallow-Water Dynamics

Shallow-Water equations describe the motion of a thin, homogeneous layer of hydrostatically-balanced, inviscid fluid with a free surface, in motion around a rotating planet (Pedlosky 1987, Holton 1992). The fluid is subject to gravitational and Coriolis forces and obeys the following equations

$`{\displaystyle \frac{𝐯}{t}}+𝐯𝐯`$ $`=`$ $`ghf𝐤\times 𝐯,`$ (1)
$`{\displaystyle \frac{h}{t}}+𝐯h`$ $`=`$ $`h𝐯,`$ (2)
where $`𝐯`$ is the horizontal velocity, $`h`$ is the thickness of the modeled layer, $`f=2\mathrm{\Omega }\mathrm{sin}\phi `$ is the Coriolis parameter, $`\mathrm{\Omega }`$ is the rotation rate of the planet, $`\phi `$ is the latitude, $`g`$ is the gravitational acceleration and $`𝐤`$ is the unit vector normal to the surface of the planet. In dimensionless form, shallow-water equations become functions only of the Rossby ($`R_o`$) and Burger ($`B_u`$) numbers:

$`R_o`$ $``$ $`{\displaystyle \frac{U}{|f|L}},`$ (3)
$`B_u`$ $``$ $`\left({\displaystyle \frac{L_D}{L}}\right)^2,L_D\sqrt{gH}/|f|,`$ (4)
where $`U`$, $`L`$ and $`H`$ are characteristic velocity, length and layer thickness scales, respectively; $`L_D`$ is the Rossby deformation radius. Note that $`|f|\mathrm{\Omega }`$ at mid-latitudes and that the planetary radius, $`R_p`$, is the relevant length scale when discussing the large-scale atmospheric circulation. The Rossby number measures the importance of rotation on the flow, while the Burger number measures the stratification of the atmosphere via the Brunt-Väisälä frequency (Holton 1992).

The atmospheric structure in bands of gaseous giant planets in our Solar System is well described as emerging from freely-evolving shallow-water turbulence on the sphere (Cho & Polvani 1996b). Turbulence in a thin atmospheric layer is quasi-2D in nature. Contrary to the forward turbulent energy cascade observed in 3D geometry, 2D turbulence is characterized by an inverse energy cascade (transfer from small to large scales) and a forward cascade of enstrophy111The flow enstrophy is defined as $`\frac{1}{2}\xi ^2`$, where $`\xi =𝐤\times 𝐯`$ is the flow vorticity. down to small scales, where it is dissipated by viscous processes. Qualitatively, the inverse cascade is associated with the growth of vortices through continuous mergers.

Turbulence in a thin atmospheric layer is also strongly constrained by the combined effect of spherical geometry and rotation (the “$`\beta `$effect”). While the force balance is everywhere the same along a latitude circle, it changes with latitude because of the dependence of the Coriolis term with $`\phi `$. This anisotropy, as measured by the parameter $`\beta =2\mathrm{\Omega }\mathrm{cos}\phi /R_p`$ (the latitudinal gradient of $`f`$), is strongest at the equator, where $`\beta `$ is maximum.. While fluid motions are free to grow to the largest available scale in the longitudinal direction, their growth is limited in the latitudinal direction by the characteristic Rhines scale, $`L_\beta =\pi \sqrt{2U/\beta }`$ (Rhines 1975). This anisotropic growth is a likely origin of the banded structure on gaseous giant planets in our Solar System; the number of bands expected for a given planet is roughly $`N_{\mathrm{band}}\pi R_p/L_\beta `$.

Cho & Polvani (1996a) presented an extensive numerical study of freely-evolving shallow-water turbulence. They explored the entire parameter space of the equations, as determined by the two dimensionless numbers $`R_o`$ and $`B_u`$. They showed that the anisotropy due to the $`\beta `$effect on a rotating sphere is necessary but not sufficient to produce a long-lasting banded structure. The Rossby deformation radius, $`L_D`$, must also be $`<R_p/3`$ for the banded structure to be stable. This small value of $`L_D`$, which acts as a limiting scale for vortex interactions, prevents the formation (via successive mergers) of large-scale structures such as circumpolar vortices.

Cho et al. (2002b) presented a generalization of the results of Cho & Polvani in the case when the planet is subject to day-side hemispheric heating, as expected for close-in EGPs. Day-side heating was prescribed in the adiabatic limit by forcing the fluid to be permanently thicker on that side, while keeping the average thickness constant. Extensive exploration of the parameter space of the forced model showed that previous shallow-water dynamics results were recovered even in the presence of this extra forcing (i.e. $`L_\beta `$ and $`L_D`$ remain the relevant scales determining the number of bands and the formation of polar vortices).

These results can be recast in terms of the values of $`R_o`$ and $`B_u`$ for the atmospheric flow. By setting $`LR_p`$, $`|f|\mathrm{\Omega }`$ and $`\beta \mathrm{\Omega }/R_p`$ (mid-latitudes), we see that the number of bands/jets expected is $`N_{\mathrm{band}}1/\sqrt{(}2R_o)`$ and that the presence of circumpolar vortices is expected for $`B_u>\mathrm{\hspace{0.33em}1}/9`$. Thus, if the values of $`R_o`$ and $`B_u`$ for other close-in EGPs can be estimated, one can get an idea of the type of large-scale atmospheric circulation expected on these planets in the stable, radiative region. In our estimates of $`R_o`$ and $`B_u`$, we will set $`U=\overline{U}`$, which is the global velocity scale of the atmospheric flow and is only known for the Solar System giants (Table 1),

An important parameter entering the definition of both $`R_o`$ and $`B_u`$ is the planetary rotation rate, $`\mathrm{\Omega }`$. While the value of $`\mathrm{\Omega }`$ is generally unknown for EGPs, a number of close-in EGPs have the advantage of being probably tidally-synchronized to their parent star, so that their rotation rate has effectively been measured via the orbital period ($`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{orb}}`$ for circular orbits). As we show below, the knowledge of $`\mathrm{\Omega }`$ for this sample of close-in EGPs restricts the range of possible values for $`R_o`$ and $`B_u`$ to a small enough region of the parameter space that their atmospheric circulation pattern can be inferred.

## 3. Sample of Close-in Extrasolar Giant Planets

Since tidal synchronization occurs faster than orbital circularization, it is possible that some close-in EGPs with substantial eccentricities are nonetheless (pseudo-)synchronized (i.e. synchronized at the periastron orbital frequency). The shallow-water results described in §2 were established only in the limit of negligible eccentricity, however. Hence, we must restrict our sample to planets with small eccentricities.222We limit the eccentricity of planets in our sample to $`e<\mathrm{\hspace{0.33em}0.1}`$, by analogy with the small eccentricities of Solar System giants to which “zero-eccentricity” shallow-water models have been applied with success (Cho & Polvani 1996b; Cho et al. 2002b). Table 1 lists all the EGPs selected for our study, plus the four Solar System giants. Parameters for the EGPs were collected from the extrasolar planet almanac333http://exoplanets.org/almanacframe.html and encyclopedia.444http://www.obspm.fr/encycl/encycl.html

Low-eccentricity EGPs were divided into two groups, based on their orbital distance to the parent star. In the first group, EGPs with semi-major axes $`a0.066`$ AU are most likely tidally-synchronized since all known EGPs with such small values of $`a`$ are also circularized. The tidal synchronization status of the more distant planets (in the second group) is less clear because several eccentric EGPs with distances of closest approach555The periastron distance, where tidal forces are the strongest, is also roughly the circular radius expected for the orbit after complete circularization. as small as $`0.05`$ AU are also known. It is thus not clear why some EGPs with periastron distances larger than this value would be tidally-circularized while others would not. We note, however, that for values of the tidal parameter $`Q`$ not too different from that of Jupiter ($`10^5`$), EGPs in this second group are also expected to be synchronized. We will assume it is indeed the case in our calculations.

Radial velocity surveys only measure $`M_p\mathrm{sin}i`$, which is a lower limit to the planet’s mass, $`M_p`$, given the unknown orbital inclination, $`i`$. For randomly oriented systems, the distribution of $`\mathrm{cos}i`$ is uniform. We adopt the value of $`M_p`$ corresponding to $`\mathrm{sin}i=0.5`$ for our fiducial estimate of $`R_o`$ and $`B_u`$, and we allow $`\mathrm{sin}i`$ to vary from $`0.1`$ to $`1`$ when estimating the range of likely values for $`R_o`$ and $`B_u`$.666Note that for low values of $`\mathrm{sin}i`$, some EGPs in Table 1 have $`M_p>13M_{\mathrm{Jup}}`$ and are thus brown dwarfs. We expect shallow-water results to be applicable even in that limit.

For a given mass, $`M_p`$, the radius of an isolated planet is estimated from the mass–radius relation for sub-stellar objects of Chabrier & Baraffe (2000), supplemented at the low mass end by a constant density law that empirically fits values for the Solar System giants. To account for the slower cooling under strong stellar irradiation, we also allow the radius to be up to $`50\%`$ larger than the value for an isolated planet, in agreement with published cooling EGP models (Burrows et al. 2000; Guillot & Showman 2002). Our results depend only weakly on the planetary radius, as long as $`R_pR_{Jup}`$ (as expected for all masses of interest). The gravitational acceleration is derived as $`g=GM_p/R_p^2`$, where $`G`$ is the gravitational constant.

For the mean layer thickness, $`H`$, we adopt the atmospheric pressure scale-height, $`H_{\mathrm{atm}}T_{\mathrm{atm}}/g`$, where $``$ is the perfect gas constant. The global radiative equilibrium temperature of the planet is $`T_{\mathrm{atm}}=T_{}(R_{}/2a)^{1/2}(1A_b)^{1/4}`$, which is a function of the parent star luminosity ($`L_{}T_{}^4R_{}^2`$), the planet’s semi-major axis $`a`$, and Bond albedo $`A_b`$. We adopt $`A_b=0.5`$ for all our numerical estimates; our results only weakly depend on the value of $`A_b`$ unless it approaches unity. The stellar luminosity is derived from the mass through the simple mass-luminosity relation $`L_{}=(M_{}/M_{})^{3.6}L_{}`$.

The last two parameters needed to determine $`R_o`$ and $`B_u`$ are the planetary rotation rate $`\mathrm{\Omega }`$ and the global kinetic energy scale $`\overline{U}`$. We assume that $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{orb}}`$ (as determined by radial velocity surveys) in all cases. We allow $`\overline{U}`$ to vary from $`50`$ m s-1, the smallest observed value for giant planets in the Solar System (Jupiter), to $`1000`$ m s-1, a rather large value for which the typical wind speeds in the atmosphere of hot, close-in EGPs approaches the sound speed. A value $`\overline{U}=400`$ m s-1 is adopted for our fiducial estimate of $`R_o`$ and $`B_u`$.

## 4. Results

Estimated values for $`R_o`$ and $`B_u`$ are given in Table 1 for Solar System giants and close-in EGPs. The values listed for close-in EGPs correspond to the range of min./max. values found given the various assumptions detailed in §3. Fiducial estimates are also reported in figure 1, where solid dots correspond to group 1 EGPs (safe tidal synchronization assumption) and open circles to group 2 EGPs (tidal synchronization assumption less safe). HD 209458 b is indicated as a star.

It is clear from figure 1 that close-in EGPs, as a group, occupy a different region of the $`R_o`$$`B_u`$ parameter space than Solar System giants. In particular, they systematically have a Burger number $`B_u>1/9`$ (even when accounting for the large range of allowed values; Table 1), which indicates that the presence of circumpolar vortices is expected in the radiative region of close-in EGPs within the framework of shallow-water dynamics. The larger values of $`R_o`$ also indicate that generally few bands/jets are expected on these planets (the uncertainty on $`\overline{U}`$ strongly affects this number; see Table 1), thus allowing the formation of larger “great spots” (which could also contribute to the variability; Cho et al. 2002b). The near alignment of all the points representing close-in EGPs in figure 1 shows that the dominant parameter determining their position in this diagram is their rotation rate ($`R_o\mathrm{\Omega }^1`$; $`B_u\mathrm{\Omega }^2`$). The small values of $`R_o`$ and $`B_u`$ for Solar System giants reflect their relatively fast rotation rates.

Although we argued in favor of variable atmospheric signatures for close-in EGPs, it is important to note that models do not yet quantitatively predict how much variability is expected. In Cho et al. (2002a,b), we emphasized that the combination of $`\overline{U}`$ (unknown) and the amplitude of day-night heating (parametrized in adiabatic simulations) determines the contrast of the thermal spots associated with circumpolar vortices. In the future, diabatic shallow-water models will allow a self-consistent determination of the day-night forcing. More sophisticated models, combined with detailed radiative transfer and chemistry descriptions (Seager & Sasselov 1998; 2000; Seager et al. 2000), will allow us to make quantitative predictions regarding the level of variability expected for various atmospheric signatures.

## Acknowledgments

Support for this work was provided by NASA through Chandra Fellowship grant PF9-10006 awarded by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073.