text:
# Spinning eggs—which end will rise?

## I Introduction

Spinning objects have historically been interesting subjects to study. The spin reversal of the rattlebackGH  (also called a celt or wobblestone) and the behavior of the tippe topGN  are typical examples. Recently, the riddle of spinning eggs was resolved by Moffatt and Shimomura.MS  When a hard-boiled egg is spun sufficiently rapidly on a table with its axis of symmetry horizontal, the axis will rise from the horizontal to the vertical. They discovered that if an axisymmetric body is spun sufficiently rapidly, a gyroscopic balance condition holds. Given this condition a constant of the motion exists for the spinning motion of an axisymmetric body. The constant, which is known as the Jellett constant,Jellett  has been found previously for symmetric tops such as the tippe top. Using these facts, they derived a first-order differential equation for $`\theta `$, the angle between the axis of symmetry and the vertical axis. For a uniform spheroid as an example they showed that the axis of symmetry indeed rises from the horizontal to the vertical.

The shape of an egg looks like a spheroid, but is not exactly so. It has thin and fat ends. Which end of the spinning egg will rise? Empirically, we know that either end can rise. But we more often see eggs spinning at the fat end with the thin end up rather than the other way round. In this paper we investigate the spinning behavior of egg-shaped axisymmetric bodies whose cross sections are described by several oval curves. We use the gyroscopic balance condition and analyze the slip velocity of the body at the point of contact as a function of $`\theta `$ and find the existence of the critical angle $`\theta _c`$ for each model curve. When the bodies are spun with the initial angle $`\theta _{\mathrm{initial}}>\theta _c`$, $`\theta `$ will increase to $`\pi `$, which means that the body will spin at the thin end. Alternatively, if $`\theta _{\mathrm{initial}}<\theta _c`$, then $`\theta `$ will decrease. For some oval curves, $`\theta `$ will decrease to 0 and the corresponding bodies will shift to the stable spinning state at the fat end. For other oval curves, a fixed point at $`\theta _f`$ is predicted, where $`0<\theta _f<\theta _c`$. In this case the bodies will spin not at the fat end but at a new stable point with $`\theta _f`$. We also explain why we observe more eggs spinning at the fat end than at the thin end.

The paper is organized as follows: To explain our notation and the geometry, we review the work of Ref. MS,  on spinning eggs in Sec. II. Then in Sec. III we introduce several models of oval curves that we will study. In Sec. IV we analyze the spinning behavior of axisymmetric bodies whose cross sections are described by these oval curves. The final section is devoted to a summary and discussion.

## II Spinning egg

We follow the geometry and notation of Ref. MS,  in their analysis of spinning eggs as much as possible. As is shown in Fig. 1, an axisymmetric body spins on a horizontal table with point of contact $`P`$. We will work in a rotating frame of reference $`OXYZ`$, where the center of mass is at the origin, $`O`$. The symmetry axis of the body, $`Oz`$, and the vertical axis, $`OZ`$, define a plane $`\mathrm{\Pi }`$, which precesses about $`OZ`$ with angular velocity $`𝛀(t)=(0,0,\mathrm{\Omega })`$. We choose the horizontal axis $`OX`$ in the plane $`\mathrm{\Pi }`$ and thus $`OY`$ is vertical to $`\mathrm{\Pi }`$ and inward. The angle of interest is $`\theta (t)`$, the angle between $`OZ`$ and $`Oz`$.

In a rotating frame of reference $`Oxyz`$, where $`Ox`$ is in the plane $`\mathrm{\Pi }`$ and perpendicular to the symmetry axis $`Oz`$ and where $`Oy`$ coincides with $`OY`$, the body spins about $`Oz`$ with the rate $`\dot{\psi }`$. Since, in the frame $`Oxyz`$, $`𝛀`$ is expressed as $`𝛀=\mathrm{\Omega }\mathrm{sin}\theta \widehat{x}+\mathrm{\Omega }\mathrm{cos}\theta \widehat{z}`$, the angular velocity of the body, $`𝝎`$, is given by $`𝝎=\mathrm{\Omega }\mathrm{sin}\theta \widehat{x}+\dot{\theta }\widehat{y}+n\widehat{z}`$. Here $`\widehat{x}`$, $`\widehat{y}`$, and $`\widehat{z}`$ are unit vectors along $`Ox`$, $`Oy`$, and $`Oz`$, respectively, $`n(t)`$ is given by $`n=\mathrm{\Omega }\mathrm{cos}\theta +\dot{\psi }`$, and the dot represents differentiation with respect to time. The $`Ox`$ and $`Oy`$ axes are not body-fixed axes but they are principal axes, so that the angular momentum, $`𝑳`$, is expressed by $`𝑳=A\mathrm{\Omega }\mathrm{sin}\theta \widehat{x}+A\dot{\theta }\widehat{y}+Cn\widehat{z}`$, where $`(A,A,C)`$ are the principal moments of inertia at $`O`$.

The coordinate system $`OXYZ`$ is obtained from the frame $`Oxyz`$ by rotating the latter about the $`Oy`$ ($`OY`$) axis through the angle $`\theta `$. Hence, in the rotating frame $`OXYZ`$, $`𝝎`$, and $`𝑳`$ have components

$`𝝎`$ $`=`$ $`((n\mathrm{\Omega }\mathrm{cos}\theta )\mathrm{sin}\theta ,\dot{\theta },\mathrm{\Omega }\mathrm{sin}^2\theta +n\mathrm{cos}\theta ),`$ (1)
$`𝑳`$ $`=`$ $`((CnA\mathrm{\Omega }\mathrm{cos}\theta )\mathrm{sin}\theta ,A\dot{\theta },A\mathrm{\Omega }\mathrm{sin}^2\theta +Cn\mathrm{cos}\theta ),`$ (2)
respectively. The evolution of $`𝑳`$ is governed by Euler’s equation,

$$\frac{𝑳}{t}+𝛀\times 𝑳=𝑿_P\times (𝑹+𝑭),$$
(3)
where $`𝑿_P`$ is the position vector of the contact point $`P`$ from $`O`$, $`𝑹`$ is the normal reaction at $`P`$, $`𝑹=(0,0,R)`$, with $`R`$ being of order $`Mg`$, the weight, and $`𝑭`$ is the frictional force at $`P`$. Because the point $`P`$ lies in the plane $`\mathrm{\Pi }`$, $`𝑿_P`$ has components $`(X_P,0,Z_p)`$, which are given by

$`Z_P`$ $`=`$ $`h(\theta )`$ (4a)
$`X_P`$ $`=`$ $`{\displaystyle \frac{dh}{d\theta }},`$ (4b)
where $`h(\theta )`$ is the height of $`O`$ above the table. We will see in Sec. IV that $`h(\theta )`$ is determined as a function of $`\theta `$, once the geometry and density distribution of the body are known.

When the frictional force is weak and $`\dot{\theta }`$ is correspondingly small, the slip velocity of the point $`P`$ is, to leading order in $`\dot{\theta }`$, expressed as $`𝑼_P=(0,V_P,0)`$, where

$$V_P=(\mathrm{\Omega }\mathrm{sin}^2\theta +n\mathrm{cos}\theta )\frac{dh}{d\theta }+(n\mathrm{\Omega }\mathrm{cos}\theta )h(\theta )\mathrm{sin}\theta .$$
(5)
Hence, the frictional force, $`𝑭`$, is to leading order, $`𝑭=(0,F,0)`$, where $`F`$ is a function of $`V_P`$ given by the law of dynamic friction between the two surfaces in contact. We assume later Coulomb friction for $`F`$.

The $`Y`$-component of Eq. (3) is expressed by

$$A\ddot{\theta }+(CnA\mathrm{\Omega }\mathrm{cos}\theta )\mathrm{\Omega }\mathrm{sin}\theta =RX_P.$$
(6)
Since the secular change of $`\theta `$ with which we are concerned is slow and thus $`|\ddot{\theta }|\mathrm{\Omega }^2`$, the first term of Eq. (6) can be neglected. Furthermore, in a situation where $`\mathrm{\Omega }^2`$ is sufficiently large so that the terms involving $`\mathrm{\Omega }`$ in Eq. (6) dominate the term $`RX_P`$, Eq. (6) is reduced, at leading order, to $`(CnA\mathrm{\Omega }\mathrm{cos}\theta )\mathrm{\Omega }\mathrm{sin}\theta =0`$. Then, for $`\mathrm{sin}\theta 0`$, we arrive at a condition

$$Cn=A\mathrm{\Omega }\mathrm{cos}\theta ,$$
(7)
which was recently discovered by Moffatt and ShimomuraMS  and was called by them as the gyroscopic balance condition. Under this condition, the Jellett constantJellett  also exists for a general axisymmetric body. With Eq. (7), the angular momentum simplifies to $`𝑳=(0,A\dot{\theta },A\mathrm{\Omega })`$, and the $`X`$\- and $`Z`$-components of Eq. (3) reduce, respectively, to

$`A\mathrm{\Omega }\dot{\theta }=FZ_P,`$ (8a)
$`A\dot{\mathrm{\Omega }}=FX_P.`$ (8b)
Equation (8), together with Eq. (4), leads to

$$𝑳𝑿_P=A\mathrm{\Omega }h=J=\text{constant},$$
(9)
where the Jellett constant $`J`$ is determined by the initial conditions. From Eqs. (4), (9), and (8a), we obtain a first-order differential equation for $`\theta `$,

$$J\dot{\theta }=Fh^2(\theta ).$$
(10)
If we assume Coulomb friction, $`F`$ is given by

$$F=\mu Mg\frac{V_P}{|V_P|},$$
(11)
and $`V_P`$ in Eq. (5), given the gyroscopic balance condition (7), is expressed as a function of $`\theta `$ as:

$$V_P=\frac{J}{Ah(\theta )}\left[\left(\mathrm{sin}^2\theta +\frac{A}{C}\mathrm{cos}^2\theta \right)\frac{dh}{d\theta }+\mathrm{sin}\theta \mathrm{cos}\theta \left(\frac{A}{C}1\right)h(\theta )\right].$$
(12)
Hence, if we know $`h(\theta )`$ from geometrical considerations , we may solve Eq. (10) and determine the time-dependence of $`\theta `$. Moffatt and ShimomuraMS  considered a uniform spheroid as an example and showed that $`\theta `$ decreases from $`\pi /2`$ to 0 for the prolate spheroid while $`\theta `$ increases from 0 to $`\pi /2`$ for the oblate one.

However, the shape of an egg is not a spheroid and has thin and fat ends. Which end of the spinning egg will rise? Empirically, we know that either end may rise and that the body spins with its axis of symmetry vertical. Which end the spinning egg chooses might seem to depend on the inclination of the initial axis of symmetry, that is, the initial value of $`\theta `$. In the following we examine several models of oval curves and determine the relation between the initial value of $`\theta `$ and the final spinning position of the egg-shaped body.

## III Models of oval curve

The shape of a three-dimensional egg can be reconstructed by rotating its two-dimensional cross section around the axis of symmetry. The cross section of an egg looks similar to an ellipse, but is not quite. It is sharper at one end than at the other. We will examine several model curves that have been proposed for the cross-section of a real egg.

Let us consider an axisymmetric body whose cross-section is described by

$$x^2=g(z),$$
(13)
with $`g(z)>0`$ for $`z_{\mathrm{min}}<z<z_{\mathrm{max}}`$ and $`g(z_{\mathrm{min}})=g(z_{\mathrm{max}})=0`$, where we choose the $`z`$ axis as the symmetry axis. If the body has uniform density, then the volume and the $`z`$ component of center of mass are given, respectively, by

$`V`$ $`=`$ $`\pi {\displaystyle _{z_{\mathrm{min}}}^{z_{\mathrm{max}}}}g(z)𝑑z,`$ (14a)
$`z_g`$ $`=`$ $`{\displaystyle \frac{\pi }{V}}{\displaystyle _{z_{\mathrm{min}}}^{z_{\mathrm{max}}}}zg(z)𝑑z.`$ (14b)
The principal moments of inertia at center of mass are expressed by

$`A`$ $`=`$ $`M{\displaystyle \frac{\pi }{V}}{\displaystyle _{z_{\mathrm{min}}}^{z_{\mathrm{max}}}}\left[{\displaystyle \frac{1}{4}}[g(z)]^2+g(z)(zz_g)^2\right]𝑑z,`$ (15a)
$`C`$ $`=`$ $`{\displaystyle \frac{M}{2}}{\displaystyle \frac{\pi }{V}}{\displaystyle _{z_{\mathrm{min}}}^{z_{\mathrm{max}}}}[g(z)]^2𝑑z.`$ (15b)
Of course, the density of a real egg is not uniform. But if the density distribution is given by $`\rho (r,z)`$ as a function of $`r`$ and $`z`$, where $`r`$ is the distance from the symmetry axis, and the cross-section is still described by Eq. (13), we can calculate $`z_g`$, $`A`$, and $`C`$.

The following are the oval curves that we will examine.

1. Cartesian oval. The curve, given by $`\sqrt{z^2+x^2}+m\sqrt{(z+a)^2+x^2}=c`$, consists of two ovals. For definiteness we set $`m=2`$. The inside oval is expressed by $`x^2=g(z)`$ with

$$g(z)=z^2\frac{8}{3}az+\frac{a^2}{9}(5\kappa ^212)\frac{4}{9}\kappa a^2\sqrt{\kappa ^23\frac{6z}{a}},$$
(16)
with $`\kappa =c/a`$. For $`\kappa =9/4`$, we find that $`g(z)`$ is defined for the interval $`\frac{17}{12}az\frac{1}{12}a`$, $`z_g=0.710a`$, and $`A/C=1.26`$ (see Fig. 2).
2. Cassini oval. This quartic curve is expressed by $`[(z+a)^2+x^2][(za)^2+x^2]=b^4`$ with $`a,b>0`$. If $`a>b`$, the curve consists of two loops, both of which look like the cross-section of a real egg with thin and fat ends. We choose the one that is expressed by $`x^2=g(z)`$ with

$$g(z)=(z^2+a^2)+a\sqrt{4z^2+\lambda ^4a^2},$$
(17)
where $`\lambda =\frac{b}{a}<1`$ for $`a\sqrt{1+\lambda ^2}za\sqrt{1\lambda ^2}`$, so that the thin end points to the positive $`z`$ axis (see Fig. 3). For $`\lambda =0.98`$, we find $`z_{\mathrm{min}}=1.40a`$, $`z_{\mathrm{max}}=0.199a`$, $`z_g=0.840a`$, and $`A/C=1.25`$.
3. Cassini oval with an air chamber. A real egg has an air chamber near the fat end. We take into account the existence of an air chamber by using the Cassini oval (17) and taking $`z_{\mathrm{min}}=\alpha a`$, with $`\sqrt{1\lambda ^2}<\alpha <\sqrt{1+\lambda ^2}`$ for the evaluation of $`V`$, $`z_g`$, $`A`$, and $`C`$. This condition means that an empty space exists for $`\sqrt{1+\lambda ^2}az\alpha a`$ (see Fig. 4). For $`\lambda =0.98`$ and $`\alpha =1.2`$, we obtain $`z_g=0.798a`$ and $`A/C=1.07`$. The position of the center of mass $`z_g`$ is closer to the thin end and the ratio $`A/C`$ is smaller compared to the curve without an air chamber.
4. Wassenaar egg curve. A rather simple equation for an oval curve was proposed recently by WassenaarWassenaar  and is given by

$$x^2=g(z)=2a\left[2z\xi a+\sqrt{4a^2+4\xi az+\xi ^2a^2}\right],$$
(18)
for $`5<\xi <6`$ and $`aza`$. For $`\xi =5.6`$ we find $`z_g=0.0714a`$ and $`A/C=1.21`$ (see Fig. 5).
5. Lemniscate of Bernoulli. The lemniscate of Bernoulli is not a candidate for oval curves and actually looks like the infinity symbol. We study it because its final spinning position might be interesting. The curve is expressed by $`[(z+a)^2+x^2][(za)^2+x^2]=a^4`$. We study the half of the curve that is given by

$$x^2=g(z)=(z^2+a^2)+a\sqrt{4z^2+a^2},$$
(19)
with $`\sqrt{2}az0`$. The lemniscate is a special case of a Cassini oval and is obtained by setting $`a=b`$ in Eq. (17). For the lemniscate, we obtain $`z_g=0.813a`$ and $`A/C=1.34`$ (see Fig. 6).

We plot in Fig. 7 the Cartesian and Cassini ovals and the Wassenaar egg curve, adjusting the parameter $`a`$ for each case so that they have the same length along the symmetry axis. We observe that the Cartesian and Cassini ovals almost overlap. Indeed the axisymmetric bodies whose cross sections are expressed by these oval curves have close values of $`A/C`$ (1.26 and 1.25 for the Cartesian and Cassini ovals, respectively). However, we will see in Sec. IV that these ovals predict different spinning behavior for the corresponding axisymmetric bodies.

## IV Which end will rise?

We obtain from Eqs. (10) and (11)

$$\dot{\theta }=\frac{\tau }{|V_P|}\stackrel{~}{V}_P,$$
(20)
with

$$\stackrel{~}{V}_P=V_P\frac{A}{J}\mathrm{and}\tau =\frac{h^2\mu Mg}{A}.$$
(21)
Equation (20) implies that the change of $`\theta `$ is governed by the sign of $`\stackrel{~}{V}_P`$. If $`\stackrel{~}{V}_P`$ is positive (negative), $`\theta `$ will increase (decrease) with time. Therefore a close examination of the behavior of $`\stackrel{~}{V}_P`$ as a function of $`\theta `$ will be important. Moffatt and ShimomuraMS  showed that for a uniform prolate spheroid, $`\stackrel{~}{V}_P`$ has the form $`\stackrel{~}{V}_P\mathrm{sin}2\theta `$ with a negative proportionality constant. Thus if the body is spun (sufficiently rapidly) on a table with the initial inclination angle $`\theta _{\mathrm{initial}}<\pi /2`$, then $`\theta `$ decreases to 0. On the other hand, if the body is spun with $`\theta _{\mathrm{initial}}>\pi /2`$, $`\theta `$ increases to $`\pi `$. Either end will rise, because both ends of the prolate spheroid look the same. The case of a real egg is different. We can easily distinguish between the thin and the fat end. We will now analyze the axisymmetric bodies whose cross-sections are expressed by the oval curves introduced in Sec. III.

We take a coordinate system in which the center of mass $`O`$ resides at the origin. In this coordinate system, the oval curves satisfy

$$x^2=f(z)=g(z+z_g),$$
(22)
where $`g(z)`$ is introduced in Eq. (13) to describe the cross-section of an axisymmetric body whose center of mass is at $`z=z_g`$. We consider the point $`P`$($`z,x=\sqrt{f(z)}`$) on the curve (see Fig. 8). The slope $`\beta `$ of the line tangent to the curve at $`P`$ is given by

$$\beta \frac{dx}{dz}=\frac{f^{}(z)}{2\sqrt{f(z)}}.$$
(23)
Draw a line from the origin which is perpendicular to the line tangent to the curve $`x=\sqrt{f(z)}`$ at $`P`$. Let the point of intersection be $`Q`$($`z_Q,x_Q`$), whose coordinates are

$`z_Q`$ $`=`$ $`{\displaystyle \frac{\beta }{\beta ^2+1}}(\beta z\sqrt{f(z)}),`$ (24a)
$`x_Q`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}z_Q.`$ (24b)
Suppose that the line $`PQ`$ is in a horizontal plane of table and $`P`$ is the point of contact. Then the line $`QO`$ defines the vertical axis $`OZ`$. The polar angle, $`\theta `$, between $`OZ`$ and $`Oz`$ is determined by,

$$\mathrm{tan}\theta =\frac{1}{\beta }=\frac{2\sqrt{f(z)}}{f^{}(z)},$$
(25)
which gives the relation between $`\theta `$ and $`z`$. The height, $`h(\theta )`$, of $`O`$ above the table is equal to the length of $`OQ`$, and we obtain

$$h(\theta )=\sqrt{z_Q^2+x_Q^2}=\frac{1}{\sqrt{\beta ^2+1}}(\sqrt{f(z)}\beta z),$$
(26)
because $`(\sqrt{f(z)}\beta z)>0`$. The squared length of $`PQ`$ corresponds to $`X_P^2`$. We choose the sign of $`X_P`$ to be the same as that of $`(zz_Q)`$ and obtain

$$X_P=\frac{1}{\sqrt{\beta ^2+1}}\left(z+\beta \sqrt{f(z)}\right).$$
(27)
If we use Eqs. (25)–(27), we confirm Eq. (4b) and find that $`V_P`$ in Eq. (12) can be rewritten as a function $`z`$ as follows:

$$V_P=\frac{J}{A}\frac{\beta ^2}{\beta ^2+1}\left[\left(\frac{1}{\beta ^2}+\frac{A}{C}\right)\frac{z+\beta \sqrt{f(z)}}{\sqrt{f(z)}\beta z}+\frac{1}{\beta }\left(\frac{A}{C}1\right)\right].$$
(28)
From Eq. (25), $`X_P`$ and $`V_P`$ can be considered as functions of $`\theta `$. As an example, we plot in Fig. 9 the graph of $`X_P`$ versus $`\theta `$ for the case of the Cartesian oval in Eq. (16). In addition to $`\theta =0`$ and $`\pi `$, $`X_P`$ vanishes at an angle $`\theta _r`$, which is obtained by solving

$$z+\beta \sqrt{f(z)}=0.$$
(29)
When the body is placed at rest on a table, its inclination angle is $`\theta _r`$ and the height $`h(\theta )`$ of center of mass $`O`$ from the table is a minimum at $`\theta _r`$. We observe from Eqs. (5) or (12) that $`V_P=0`$ at $`\theta =0`$ and $`\pi `$, because $`\mathrm{sin}\theta =0`$ and $`\frac{dh}{d\theta }(=X_p)=0`$ at these points. Moreover, $`V_P`$ vanishes at other angles, which are given by solving

$$\frac{A}{C}+\frac{z}{\beta \sqrt{f(z)}}=0.$$
(30)
When $`A=C`$, Eqs. (29) and (30) become equivalent, which means that $`V_P`$ and $`X_P`$ vanish at the same inclination angle.

We next examine the graph of $`\stackrel{~}{V}_P`$ as a function of $`\theta `$ for the oval curves introduced in Sec. III.

(i) A Cartesian oval. Figure 10 shows that $`\stackrel{~}{V}_P`$ crosses the line $`\stackrel{~}{V}_P=0`$ at an angle $`\theta _c`$ and $`\stackrel{~}{V}_P>0`$ for $`\theta _c<\theta <\pi `$ but is negative for $`0<\theta <\theta _c`$. So, the angle $`\theta _c`$ is a critical point. If the body is spun on a table with the initial angle $`\theta _{\mathrm{initial}}>\theta _c`$, then $`\theta `$ will increase to $`\pi `$, which means that the body will eventually spin at the thin end. For $`\theta _{\mathrm{initial}}<\theta _c`$, we will see that the body spins at the fat end. That is, depending on the initial value $`\theta _{\mathrm{initial}}`$, the body will spin at the thin or the fat end. Both ends are stable points. We have $`A>C`$ for Cartesian ovals, which leads to $`\theta _c>\theta _r`$. Numerically we obtain from Eqs. (30) and (29) that $`\theta _c=1.86`$ and $`\theta _r=1.72`$. Recall that $`\theta _r`$ is the inclination angle when the body is placed at rest. If we spin the body involuntarily, the initial angle $`\theta _{\mathrm{initial}}`$ tends to be near $`\theta _r`$. Because $`0<\theta _r<\theta _c`$, it is likely that the body is spun with $`\theta _{\mathrm{initial}}`$ between 0 and $`\theta _c`$, and thus it will shift to the stable spinning state at the fat end. Empirically, we more often observe eggs spinning at the fat end rather than at the thin end. The expected behavior of the axisymmetric body expressed by the Cartesian oval in Eq. (16) well explains the observed features of the spinning egg.

(ii) A Cassini oval. The second example of oval curves presents an interesting situation. We see from Fig. 11 that $`\stackrel{~}{V}_P`$ crosses the line $`\stackrel{~}{V}_P=0`$ at $`\theta _f`$ and $`\theta _c`$ ($`\theta _f<\theta _c`$), and that $`\stackrel{~}{V}_P`$ is negative for $`\theta _f<\theta <\theta _c`$ and otherwise positive. Numerically we obtain $`\theta _f=0.45`$ and $`\theta _c=1.92`$ (and $`\theta _r=1.77`$). Thus the graph of $`\stackrel{~}{V}_P`$ in Fig. 11 implies that for the Cassini oval (17) the thin end ($`\theta =\pi `$) is a stable point, but the fat end is not. When the body is spun with the initial value of $`\theta `$ anywhere between 0 and $`\theta _c`$, $`\theta `$ will approach the fixed point $`\theta _f`$. In other words, the body will spin not at the fat end, but at the point with the inclination angle $`\theta _f`$. It is interesting to note that the curves of the Cartesian (16) and Cassini (17) ovals almost overlap each other when they are adjusted to have the same length along the symmetry axis (see Fig. 7). But they predict different behaviors for $`\stackrel{~}{V}_P`$ and thus different spinning behaviors for the corresponding axisymmetric bodies. Does a hard boiled egg show the behavior predicted by this Cassini oval? We are almost certain that we have never seen such a behavior.

(iii) A Cassini oval with an air chamber. Because an egg has an air chamber near the fat end, we study the case of a Cassini oval with an air chamber. The existence of an air chamber moves the position of center of mass $`z_g`$ toward the thin end, from $`z_g=0.840a`$ to $`z_g=0.798a`$, and reduces the ratio $`A/C`$ from 1.25 to 1.07. Consequently, the fixed point at $`\theta _f`$, which is present for Eq. (17), disappears. The graph of $`\stackrel{~}{V}_P`$ in Fig. 12 shows that it crosses the line $`\stackrel{~}{V}_P=0`$ only once at $`\theta _c=1.97`$. The inclination angle at rest becomes $`\theta _r=1.93`$, which is very close to the value of $`\theta _c`$ because $`A/C1`$. The axisymmetric body described by a Cassini oval with an air chamber also reproduces the features of the spinning egg.

(iv) Wassenaar egg curve. Figure 13 shows that this curve has a fixed point at $`\theta _f=0.93`$. When $`\theta _{\mathrm{initial}}`$ is between 0 and $`\theta _c=1.85`$, $`\theta `$ will move to $`\theta _f`$, another stable point in addition to the one at the thin end ($`\theta =\pi `$). In addition, the graph of $`X_P`$ for this oval curve vanishes at two other points. One is at $`\theta _r=1.78`$, a position at rest, and one at $`\theta =0.58`$, an unstable point.

(v) Lemniscate of Bernoulli. For the lemniscate of Bernoulli, we find $`lim_{zz_{\mathrm{max}}}(1/\beta )=1`$, so that Eq. (25) tells us that the allowed region of $`\theta `$ is between 0 and $`3\pi /4`$. Figure 14 shows that $`\stackrel{~}{V}_P`$ vanishes at the fat end ($`\theta =0`$), at the fixed point ($`\theta _f=0.74`$), and at the critical point ($`\theta _c=1.96`$). The position of the body at rest is at $`\theta _r=1.81`$, and its spinning state has two stable points at $`\theta =\theta _f`$ and $`3\pi /4`$.

## V Summary and Discussion

We have examined the spinning behavior of axisymmetric bodies whose cross sections are described by several model curves, including a Cartesian oval, Cassini ovals with and without an air chamber, and the Wassenaar egg curve. These results together with the lemniscate of Bernoulli are summarized in Table 1. For each oval curve we used the gyroscopic balance condition (7) and found the predicted slip velocity $`V_P`$ of the contact point as a function of the inclination angle $`\theta `$ and the existence of the critical angle $`\theta _c`$. When the body is spun on a table with the initial angle $`\theta _{\mathrm{initial}}>\theta _c`$, $`\theta `$ will increase to $`\pi `$, which means that the body will spin at the thin end. If $`\theta _{\mathrm{initial}}<\theta _c`$, then $`\theta `$ will decrease. For the Cartesian oval and Cassini oval with an air chamber, $`\theta `$ will reduce to 0 and the corresponding bodies will spin at the fat end. Moreover, when the bodies are spun without intention, we expect to see their spinning states at the fat end more often than at the thin end because the inclination angle $`\theta _r`$ at rest is smaller than $`\theta _c`$. This behavior is consistent with the features of a spinning egg.

On the other hand, the Cassini oval and Wassenaar egg curves predict the existence of the fixed point at $`\theta _f`$, where $`0<\theta _f<\theta _c`$. Then the fat end ($`\theta =0`$) is no longer a stable point. If the corresponding bodies are spun with $`\theta _{\mathrm{initial}}<\theta _c`$, $`\theta `$ moves to $`\theta _f`$ and not to 0, and the bodies will spin at a new stable point at $`\theta _f`$. The lemniscate of Bernoulli is not an oval curve, but the body described by this curve also has a fixed point.

It would be very interesting to make axisymmetric bodies whose cross sections are described by the Cassini oval (17), Wassenaar egg curve (18) and the lemniscate of Bernoulli (19) and determine if those bodies will spin at a new stable point that is different from the fat end.

The inclusion of an air chamber at the fat end tends to diminish the appearance of the fixed point at $`\theta _f`$. The case of a Cassini oval is an example. For the Wassenaar egg curve, we need a rather large air chamber; an empty space for $`az0.5a`$ in Eq. (18) is necessary for the disappearance of $`\theta _f`$. The fixed point of the lemniscate vanishes if we take $`z_{\mathrm{min}}=1.33a`$ for Eq. (19). Moreover, the inclusion of an air chamber at the fat end moves the position of the center of mass $`z_g`$ toward the thin end and makes the ratio $`A/C`$ smaller. Consequently, the values of the critical angle $`\theta _c`$ and the inclination angle at rest $`\theta _r`$ move toward $`\pi `$.

We have assumed Coulomb’s law for the friction $`F`$ (see Eq. (11)). If we instead assume a viscous friction law,

$$F=\stackrel{~}{\mu }MgV_P,$$
(31)
all our conclusions remain unchanged, in particular, the positions of the critical points and fixed points. Only the transition time from the unstable to the stable state will be modified. As an example, the transition time from the angle $`\theta _c`$ to $`\pi `$ is numerically calculated to be $`t(\theta _c\pi )=(J/\mu Mga^2)\chi `$ with a numerical factor $`\chi 𝒪(1)`$ for Coulomb friction, and $`t(\theta _c\pi )=(1/\stackrel{~}{\mu }g)\chi `$ with $`\chi 𝒪(10)`$ for viscous friction .

###### Acknowledgements.

The author thanks Sinsuke Watanabe and Tsuneo Uematsu for valuable information on the spinning egg and the tippe top and for critical comments. He is indebted to Tsuneo Uematsu for reading the manuscript and important suggestions. Thanks also are due to Shingo Ishiwata, Tetsuji Kuramoto, and Yoshihiro Shimazu for useful discussions, and Takahiro Ueda for assistance with figures. This paper is dedicated to the memory of Sozaburo Sasaki.

## Figures


format: text

text:
# Spin Blockade in Capacitively Coupled Quantum Dots

## Abstract

We present transport measurements on a lateral double dot produced by combining local anodic oxidation and electron beam lithography. We investigate the tunability of our device and demonstrate, that we can switch between capacitive and tunnel coupling. In the regime of capacitive coupling we observe the phenomenon of spin blockade in a magnetic field and analyze the influence of capacitive interdot coupling on this effect.

In recent decades great progress has been made in the development and measurement of quantum dot devices Kouwenhoven et al. (1997). Many single electron transistors have been fabricated and during the last years coupled quantum dots (e.g. van der Wiel et al. (2003); Holleitner et al. (2002); Livermore et al. (1996); Hofmann et al. (1995); Molenkamp et al. (1995)) came more and more into the focus of research, as they are proposed as crucial parts for quantum computers Loss and DiVincenzo (1998). Since preparation and detection of electron spin is another essential element in future quantum information applications, the discovery of spin blockade in quantum dot devices was a great step forward Ciorga et al. (2000, 2002). It has been observed in single quantum dots defined by metallic gates. Coupled quantum dots are also of great interest in the regime of spin blockade. Only recently spin resolved measurements on molecular states of two dots tunnel coupled in series have been made using spin blockade Pioro-Ladriere et al. (2003).

In this letter we present our results for spin blockade in a parallel double quantum dot based on local anodic oxidation (LAO). Due to the complexity of the device we first have to characterize our system in terms of the interdot coupling, which depends on top gate voltage and magnetic field. We find, that the relevant regime, where we observe spin blockade, features capacitive coupling only. We show our measurements for spin blockade and analyze the influence of capacitive interdot coupling on this effect.

Our quantum dot device is based on a GaAs/AlGaAs heterostructure with a two-dimensional electron system (2DES) 34 nm below the surface. The sheet density is $`n=4.310^{15}`$ m-2. We apply two different nanolithographic techniques to the sample surface to define two adjacent quantum dots connected in parallel to common source and drain contacts. We use an atomic force microscope (AFM) to write the basic double dot structure by local anodic oxidation. We complete the structure with a metallic gate patterned with electron beam lithography (e-beam) to add the function of controlled tunability of the interdot coupling.

An AFM image of our device is given in Fig. 1(a). Those parts, that have been fabricated by LAO, are colored blue. They represent insulating lines, that form a complete double dot device including the two dots A and B, side gates G1 and G2 and source (S) and drain (D) leads. The two dots are connected via 80 to 100 nm wide tunneling barriers to the leads and via a central barrier to each other. An additional metallic top gate (golden) can be biased to tune the height of the central barrier and thus to tune the interdot coupling. A detailed description of a similar device and the fabrication technique can be found in Rogge et al. (2003).

Our double dot device is characterized in transport measurements in a 3He/4He dilution refrigerator at a base temperature of 70 mK. We apply standard Lock In technique to measure the differential conductance $`G`$ through the double dot system. To analyze the features for spin blockade, which will be shown later, it is first necessary to understand the interaction of both dots in the relevant parameter space. Therefore we start investigating the interdot coupling with the top gate grounded at $`B`$=0 T. Figure 1(c) shows $`G`$ under these conditions as a function of both side gate voltages $`V_{G1}`$ and $`V_{G2}`$. The dark regions denote high values for $`G`$ and represent transport over dot A. No lines corresponding to transport over dot B are visible meaning that B is not connected to both leads. Nevertheless the equidistant splittings interrupting the features of dot A indicate the presence of dot B, which is connected to the source lead while the drain lead remains closed for all the parameters used here, as we know from nonlinear measurements (not shown here). This is due to the lithographic gap width of only 70 nm for the barrier between dot B and drain compared to 80 to 100 nm for the other barriers. Thus we can complete the observed features by dotted lines to a hexagonal pattern typical for two coupled quantum dots. Crossing a dark line changes the electron number on dot A ($`N_A`$), crossing a dotted line changes the charge on dot B ($`N_B`$). Due to the sharp corners of the hexagons and the fact that the only features visible are those of dot A, the interdot coupling in this regime is purely capacitive.

Figure 1(d) shows a similar plot but with the top gate voltage set to 80 mV. Due to the reduced height of the central barrier at $`V_{TG}`$=80 mV the interdot coupling has increased leading to a very wide splitting. In contrast to Fig. 1(c) the edges of the hexagonal pattern are clearly rounded and there are finite conductance values visible for transitions on dot B. Both effects demonstrate the molecular like character of the system. Instead of atomic transitions the visible features refer to molecular binding states indicating tunnel coupling Blick et al. (1996, 1998).

A more detailed view can be gained by calculating the fractional peak splitting $`F=2\mathrm{\Delta }V_S/\mathrm{\Delta }V_P`$ (see Fig. 1(c)), which is 0 for two separate dots and 1 for two totally merged dots Golden and Halperin (1996). Figure 1(d) shows $`F`$ as a function of top gate voltage for $`B`$=0 T (open circles). For voltages below 30 mV the fractional peak splitting is unaffected by $`V_{TG}`$ and remains at values below 0.2. This is due to capacitive coupling, that depends on the LAO geometry of both dots, which is almost not affected by top gate voltage. At 30 to 40 mV tunnel coupling sets in and molecular features appear. With increasing top gate voltage the fractional peak splitting can be tuned to almost 0.8 corresponding to very strong tunnel coupling.

The interdot coupling is not only changed by top gate voltage but also by a magnetic field applied perpendicular to the 2DES. At $`B`$=3.3 T the fractional peak splitting is again decreased to a value below 0.2 and remains unaffected by top gate voltage (squares in Fig. 1(d)). The magnetic field destroys the tunnel coupling and moves the system back to the regime of capacitive coupling even with $`V_{TG}`$=80 mV.

After having characterized our double dot we can now focus on the features of spin blockade, that appear around 3.3 T in the regime of capacitive coupling.

Figure 2(a) shows $`G`$ as a function of $`V_{G1}`$ and $`B`$. Two Coulomb peaks are visible denoting transitions on dot A. For both peaks oscillations in peak position are obvious as well as alternating amplitudes for the peak maxima (see Fig. 2(b)). The amplitude is changed by magnetic field and by electron number $`N_A`$. Therefore this effect looks similar to the effect of spin blockade that has been observed in several single dot devices Ciorga et al. (2000, 2002).

Spin blockade is explained as follows. As source and drain consist of a 2DES, they develop edge channels in a magnetic field. If the potential is flat enough these channels are separated in space leaving a spin down channel of the lowest Landau level near to the dot and a spin up channel further away. Thus the overlap of electronic wave functions of leads and dot is spin dependent. Spin up electrons couple weaker to the dot than spin down electrons. Therefore transport with spin up electrons is suppressed. A magnetic field changes the available spin state for electrons tunneling into the dot McEuen et al. (1992), which becomes obvious in alternating peak positions and peak amplitudes of Coulomb blockade peaks. LAO devices are known to have steeper potentials than split gate devices Fuhrer et al. (2001); Held et al. (1999). Thus the spatial separation of edge channels in the leads is believed to be much smaller in LAO devices. Spin blockade should be very improbable or at least suppressed in comparison to split-gate devices and has actually not been measured in LAO devices so far.

Nevertheless, despite the LAO nature we observe these features here very clearly, and although this is a coupled dot system, we assume, that this observation is due to spin blockade in dot A only, since we are in the regime of capacitive coupling. Nevertheless the influence of dot B, which is still present, must be investigated. Therefore we measured $`G`$ as a function of magnetic field and top gate voltage. In contrast to side gate 1, that couples mainly to dot A, the top gate has a stronger influence on dot B. Thus we can expect to see features from both dots at once.

The measurement is shown in Fig. 3(a). Again we observe oscillating peak positions in the Coulomb peaks of dot A. But a detailed look reveals, that these peaks are split by local minima in differential conductance. These minima are colored white in the highlighted section and show a pattern of more or less horizontal lines with regular distance in $`V_{TG}`$. These lines refer to transitions on dot B and correspond to the splitting of triple points in Fig. 1(b). As the Coulomb peaks of dot A are split, parts of each peak are above a splitting and others are below (red lines in the section). Therefore we can compare the peak amplitudes of each Coulomb peak below and above a splitting and can directly investigate the effect of capacitive dot-dot coupling on the peak amplitude. This is done in Fig. 3(b). For each peak visible in Fig. 3(a) the peak amplitude is plotted as a function of magnetic field. For all peaks both amplitudes below a splitting (solid line) and above a splitting (dotted line) fit almost perfectly together and the typical oscillations of high and low amplitude become visible. Thus the influence of capacitive interdot coupling is identified as a pure peak splitting effect without any impact on the peak amplitude. In consequence we find our assumption verified, that the oscillations in peak position and amplitude ascribe to a single dot spin blockade of dot A. As a result this spin blockade is not destroyed by capacitive interdot coupling, it is not even affected at all.

In summary we have investigated spin blockade in a tunable double dot device fabricated by LAO. We have characterized the system regarding the tunability of the interdot coupling as a function of top gate voltage and magnetic field. In the regime of capacitive coupling we observed single dot spin blockade. As far as we know this is the first time that this has been measured in devices fabricated by LAO. We have shown that this spin blockade is not destroyed or at least affected by transitions on the capacitively coupled dot. Therefore it is a very useful tool for spin detection even in LAO based quantum dot devices.

We thank M. Bichler, G. Abstreiter and W. Wegscheider for the heterostructure and D. Pfannkuche for helpful discussions. We thank F. Hohls for careful reading of the manuscript. This work has been supported by BMBF.


text:
# THE BOLTZMANN ENTROPY FOR DENSE FLUIDS NOT IN LOCAL EQUILIBRIUM

## I INTRODUCTION

In 1854 Clausius introduced the notion of the entropy of a macroscopic system (defined up to additive constants) in a state of thermal equilibrium Jaynes1 . Not long afterwards Boltzmann gave a microscopic definition of the entropy $`S\left(M\right)`$ of a general macroscopic system in a macrostate defined by values of macrovariables $`M`$: $`M=M_X`$ represents suitable “coarse-grained” functions of the system’s microstate, given by a point $`X`$ in the $`2Nd`$-dimensional phase space $`\mathrm{\Gamma }`$ for a $`d`$-dimensional classical system containing $`N`$ particles. $`S\left(M\right)`$ is (up to constants) equal to the log of the volume $`\left|\mathrm{\Gamma }_M\right|`$ of the phase-space region $`\mathrm{\Gamma }_M`$, defined by the macrostate $`M`$, i.e. containing all phase points (or microstates) giving rise to this macrostate Jaynes2 . $`S\left(M\right)`$ agrees with the Clausius entropy for systems in equilibrium when $`M=M_{eq}`$ is just the total energy $`E`$, of $`N`$ particles in a volume $`V`$. The fact that $`\left|\mathrm{\Gamma }_{M_{eq}}\right|`$ is exponential in the number of molecules in the system explains the origin of the second law in the microscopic dynamics: when a constraint is lifted in a system in equilibrium, thereby affording access to a new equilibrium macrostate $`M_{eq}^{}`$ of larger entropy per molecule than the original $`M_{eq}`$, the overwhelming majority of phase points in $`\mathrm{\Gamma }_{M_{eq}}`$ will find themselves eventually in $`\mathrm{\Gamma }_{M_{eq}^{}}`$, since for a macroscopic system $`\left|\mathrm{\Gamma }_{M_{eq}^{}}\right|`$ is enormously larger than $`\left|\mathrm{\Gamma }_{M_{eq}}\right|`$ as well as the volume of the union of all nonequilibrium macrostates without the constraint Jaynes2  Boltz .

Boltzmann’s interpretation of entropy naturally extends the second law to nonequilibrium macroscopic systems: a (significant) violation of the second law will not occur provided the microstate $`X`$ of a physical system prepared in, or evolved into, the macrostate $`M`$ is typical (or at least not atypical) of points in $`\mathrm{\Gamma }_M`$ as far as the future evolution of $`M`$ is concerned. Of course to rigorously prove that $`X`$ should be (and remain) typical in this sense is difficult. It was pointed out, however, in Resib  that a sufficient condition for this is that the evolution of $`M_t`$ be given by an autonomous deterministic law, i.e. the value of $`M`$ at any time $`t+\tau `$ is determined by its value at time $`t`$, $`M_{t+\tau }=\mathrm{\Phi }_\tau \left(M_t\right)`$, for $`\tau 0`$. Such a law implies that, with rare exception, $`\varphi _\tau \mathrm{\Gamma }_{M_t}\mathrm{\Gamma }_{M_{t+\tau }}`$, i.e. the overwhelming majority of phase-points (almost all in suitable limits) $`X\mathrm{\Gamma }_{M_t}`$ will, when evolved according to the microscopic evolution law $`\varphi _\tau `$, be found in $`\mathrm{\Gamma }_{M_{t+\tau }}`$ for $`\tau 0`$. It follows then from Liouville’s theorem that (with insignificant error) $`S\left(M_{t+\tau }\right)S\left(M_t\right)`$.

The most common example of a deterministic macro-evolution occurs for simple systems in local thermal equilibrium (LTE), described by the macrovariables $`M`$ representing the locally conserved (particle, momentum and energy) densities. $`M_t`$ then satisfies deterministic hydrodynamic-type equations, e.g. the Euler or Navier-Stokes equations, and $`S\left(M_t\right)`$, which is given by an integral over the volume $`V`$ of the equilibrium entropy density, see Alder , is then indeed monotone nondecreasing in time.

Such a deterministic evolution is however not necessary for the monotonicity of $`S\left(M_t\right)`$. There is an enormous disparity between the (small) number of possible macrostates (which are always defined with some macroscopic tolerance in terms of a relatively small number of macrovariables) and the large number of microstates corresponding to the possible values of the large number of microvariables (counted say in terms of phase-space cells of volume $`\overline{h}^{dN}`$). Even when $`M_0`$ does not determine $`M_t`$, $`S\left(M_{X_t}\right)`$ should be nondecreasing to leading order in the size of the system for the overwhelming majority of initial microstates $`X_0`$; see also Ehrenfest .

The choice of macrovariables $`M`$ and the corresponding computation of $`S\left(M\right)`$ suitable for describing systems not in LTE is a daunting task, especially for complex systems such as polymeric fluids, metals with memory, etc. The first and still paradigmatic step in that direction was taken by Boltzmann himself when he computed $`S\left(f\right)`$ for the macrovariables $`f=\left\{f_X(\underset{¯}{x},\underset{¯}{v})\right\}`$ corresponding to the coarse-grained (empirical) density in the $`2d`$ dimensional $`\mu `$-space for a macroscopic system in a microstate $`X`$ Jaynes2 ,Goldstein . He found that (up to constants)

$$S\left(f\right)=S_{gas}\left(f\right)k_V𝑑\underset{¯}{x}_{R^d}𝑑\underset{¯}{v}f(\underset{¯}{x},\underset{¯}{v})\mathrm{ln}f(\underset{¯}{x},\underset{¯}{v})$$
(1)
We put the subscript “gas” on $`S_{gas}\left(f\right)`$ to emphasize that $`f(\underset{¯}{x},\underset{¯}{v})`$ can be expected to suffice for the adequate specification of the macrostate away from LTE only for a dilute gas, where interactions between particles make a negligible contribution to the energy of the system. In fact it was for such a dilute gas that Boltzmann derived a deterministic evolution equation for $`f`$ and proved that the corresponding $`f_t`$ satisfies the $``$-theorem, $`\frac{d}{dt}S_{\mathrm{gas}}\left(f_t\right)0`$, Jaynes2 , Boltz , Goldstein . For dense fluids, specification of $`f(\underset{¯}{x},\underset{¯}{v})`$ is compatible with many different total energies (including infinite ones for hard core interparticle potentials). A simple analysis then shows that the phase points $`X`$ of a system with specified energy $`E`$ which is below the maximal energy compatible with $`f`$ will correspond to an exceedingly small minority of the phase points in $`\mathrm{\Gamma }_f`$, i.e. will be atypical of points in $`\mathrm{\Gamma }_f`$. There is then no reason to expect for such systems that $`S_{gas}\left(f_{X_t}\right)`$ will increase as the system evolves in time according to its energy conserving Hamiltonian dynamics. It is in fact easy to set up in dense fluids initial macrostates $`f_0`$ such that for $`X_0\mathrm{\Gamma }_{M_0}`$ $`S_{gas}\left(f_{X_t}\right)`$ will typically decrease in time when the fluid goes to equilibrium Lebowitz1 .

It was argued in Resib  that if one includes in $`M`$, in addition to $`f`$, also the total energy $`E`$, then the entropy $`S(f_t,E)`$ should be an increasing function of time, i.e. $`S(f_t,E)`$ should satisfy an $``$-theorem for general systems, including dense fluids. It was also noted there that the quantity shown by Resibois Rapaport  to satisfy an $``$-theorem for $`f_t`$ evolving via the modified Enskog equation (expected to be accurate for moderately dense hard sphere gases) is in fact the Boltzmann entropy $`S(f_t,E)`$ for a system of hard spheres.

In this work we use molecular dynamics to investigate the time evolution of $`S(f_t,E)`$ for dense fluids interacting with Lennard-Jones and other types of pair potentials. We consider in particular situations, such as those in Lebowitz1 , where $`S_{gas}\left(f_t\right)`$, defined in (1), is expected to decrease. Our simulations, which give a monotone increase of $`S(f_t,E)`$ when the number of particles in the system is large, support the hypothesis that the time evolution of a typical microstate in $`\mathrm{\Gamma }_{f,E}`$ is indeed such that $`S(f_t,E)`$ satisfies an $``$-theorem. We also find evidence that $`f_t`$ itself evolves in a deterministic way, with different microstates with the same $`f_0`$ give rise to the same $`f_t`$, although no equation yielding this evolution is at present known (at least to us) for general dense fluids. This suggests looking for an autonomous equation for $`f_t`$ (at a given $`E`$). This is exactly what is done in the heuristic derivations of the Boltzmann and Enskog equations Lennard  and is discussed extensively in the literature for various other systems, see 11 . The validity of the Boltzmann equation for dilute gases, i.e. for typical $`X_0\mathrm{\Gamma }_{f_0}`$ in the Boltzmann-Grad limit, was justified rigorously at least for short times, by Lanford Lanford .

Formalism: We consider a system of $`N`$ particles with unit mass in a box $`V`$. The microstate is specified by $`X=(\underset{¯}{x}_1,\underset{¯}{v}_1,\mathrm{},\underset{¯}{x}_N,\underset{¯}{v}_N)`$ and the dynamics is given by the Hamiltonian

$$H\left(X\right)=\frac{1}{2}\underset{i=1}{\overset{N}{}}\underset{¯}{v}_i^2+\frac{1}{2}\underset{ij}{}\varphi \left(\underset{¯}{x}_i\underset{¯}{x}_j\right)$$
(2)
It follows from the structure of the classical phase space $`\mathrm{\Gamma }`$ that $`S(f,E)`$ can be written, c.f. Resib , as

$$S(f,E)=S^{\left(m\right)}\left(f\right)+S^{\left(c\right)}(n,\mathrm{\Phi }_{total}).$$
(3)
Here $`S^{\left(m\right)}`$ ($`S^{\left(c\right)}`$) is the log of the momentum (configuration) space volume corresponding to the macro-state $`M=(f,E)`$,

$$S^{\left(m\right)}\left(f\right)=S_{gas}\left(f\right)+k_V𝑑\underset{¯}{x}n\left(\underset{¯}{x}\right)\mathrm{log}n\left(\underset{¯}{x}\right)$$
(4)
with $`n\left(\underset{¯}{x}\right)=_{R^d}𝑑\underset{¯}{v}f(\underset{¯}{x},\underset{¯}{v})`$, the spatial density, and

$$S^{\left(c\right)}(n,\mathrm{\Phi }_{tot})=\underset{\mathrm{\Phi }}{sup}_V𝑑\underset{¯}{x}s^{\left(c\right)}(n\left(\underset{¯}{x}\right),\mathrm{\Phi }\left(\underset{¯}{x}\right))$$
(5)
where $`s^{\left(c\right)}(n,\mathrm{\Phi })`$ is the configurational entropy density of an equilibrium system with Hamiltonian (2) having particle density $`n`$ and potential energy density $`\mathrm{\Phi }`$. The $`sup`$ in eq. (5) is taken over all $`\mathrm{\Phi }\left(\underset{¯}{x}\right)`$ such that

$$_V𝑑\underset{¯}{x}\mathrm{\Phi }\left(\underset{¯}{x}\right)=\mathrm{\Phi }_{total}=E_V𝑑\underset{¯}{x}_{R^d}𝑑\underset{¯}{v}f_t(\underset{¯}{x},\underset{¯}{v})\frac{1}{2}\underset{¯}{v}^2$$
(6)
Restricting ourselves to spatially uniform systems, $`n=N/V`$, $`f=nh\left(\underset{¯}{v}\right)`$, $`_{R^d}𝑑\underset{¯}{v}\underset{¯}{v}h\left(\underset{¯}{v}\right)=0`$, $`\mathrm{\Phi }=\mathrm{\Phi }_{total}/V`$, we find (see eqs. (39)-(41) in Goldstein ) that

$$\frac{d}{dt}S^{\left(c\right)}(n,\mathrm{\Phi }_{\mathrm{total}})=\frac{1}{T_\mathrm{\Phi }}\frac{d\mathrm{\Phi }_{total}}{dt}$$
(7)
where $`T_\mathrm{\Phi }`$ is the inverse of $`\mathrm{\Phi }_{eq}\left(T\right)`$, the potential energy density of the equilibrium system with Hamiltonian (2) at density $`n`$.

Simulations: To check whether $`S(f_t,E)`$ as expressed in eqs. (3) to (7) satisfies an $``$-theorem for dense fluids we have carried out simulations on a two-dimensional system with density $`n=0.5`$ in a periodic box interacting with a cut-off Lennard-Jones potential, $`\varphi \left(r\right)=\overline{\varphi }\left(r\right)\overline{\varphi }\left(r_c\right)\left(rr_c\right)d\overline{\varphi }\left(r\right)/dr|_{r=r_c}`$ with $`\overline{\varphi }\left(r\right)=4\left[r^{12}r^6\right]`$ for $`rr_c=2.5`$ and $`\varphi \left(r\right)=0`$ otherwise. To obtain $`T_\mathrm{\Phi }`$ in (7) we first computed $`\mathrm{\Phi }_{eq}\left(T\right)`$, for such a system. For these simulations, as well as for those described below, we used the Verlet algorithm with time mesh $`10^4`$.

To carry out the time dependent simulations, we first let the system reach an equilibrium state. We then multiplied the speed of each particle with an appropriate factor to obtain a state with too high a kinetic energy, and let the system evolve freely to its new equilibrium state. This means that the initial $`f`$ was a Maxwellian with too high a temperature for the total energy—the case considered in Lebowitz1 . During the consequent evolution the velocity distribution stayed isotropic and we computed $`S^{\left(m\right)}`$, as well as the kinetic and potential energy densities, in time intervals of size $`0.001`$. To compute $`S^{\left(m\right)}`$ we first find the speed distribution by counting the number of particles having their speeds in each of the $`50`$ equal intervals into which we divide the segment $`[0,v_{max}]`$ where $`v_{max}`$ is the maximum speed of any particle at time $`t`$. We then get the configurational entropy $`S^{\left(c\right)}`$ by taking a numerical derivative of $`\mathrm{\Phi }`$ with respect to time and then carrying out a numerical integration of the right side of (7). In Figure 1 we show the time behavior of $`S^{\left(m\right)}`$, $`\mathrm{\Phi }_{\mathrm{total}}`$ and $`S`$ for such a system.

We have carried out similar calculations for the truncated repulsive potential $`\varphi \left(r\right)=\overline{\varphi }\left(r\right)\overline{\varphi }\left(r_c\right)\left(rr_c\right)d\overline{\varphi }\left(r\right)/dr|_{r=r_c}`$ with $`\overline{\varphi }\left(r\right)=r^6`$ for $`rr_c=2.5`$ and $`\varphi \left(r\right)=0`$ otherwise. We show the results for different values of $`N`$ in Fig. 2. The fluctuations for small $`N`$ are clearly visible: their magnitude, once the system has reached equilibrium, appear to scale as $`N^{1/2}`$. There are also finite-size corrections to the equilibrium time averages consistent with those expected from using a micro-canonical ensemble.

We also investigated for this system whether $`f_t`$ evolves deterministically by comparing $`f_t`$ for different initial microstates, all having the same $`f_0`$ and $`E`$. The results are shown in Figs. 3 and 4. The initial $`f_0`$ in Fig. 3 is close to a Maxwellian and the subsequent $`f_t`$ are also close to a Maxwellian with time dependent temperatures. In Fig. 4 the initial $`f_0`$ is one in which all the particles have the same speed. The evolution of $`f_t`$ towards a Maxwellian is clearly visible.

Binary Mixtures: We also considered the case of a binary system of hard point particles with alternating masses $`m_1`$ and $`m_2`$ 2mass . We used as our macrovariable the total energy $`E`$ plus the empirical density $`f\left(v\right)dv`$ of particles with velocities in some interval $`dv`$ (independent of the species) with uniform positional densities $`n_1`$ and $`n_2`$ in a box of length $`L`$. The entropy $`S_{\left(2\right)}(f,E)`$ of this system can be written as a sum:

$$S_{\left(2\right)}(f,E)=S_{\mathrm{gas}}\left(f_1\right)+S_{\mathrm{gas}}\left(f_2\right),E=E_1+E_2$$

The maximum of $`S_{\left(2\right)}(f,E)`$ is obtained for $`f`$ a sum of Maxwellians $`f_1`$ and $`f_2`$ with the same temperature $`T`$ determined by the total energy $`E=\left(n_1+n_2\right)\frac{1}{2}kTL`$ .

Starting with an initial microstate for which $`S_{\mathrm{gas}}\left(f_1\right)>Ls_{\mathrm{eq}}(n_1,T_1)`$, $`S_{\mathrm{gas}}\left(f_2\right)<Ls_{\mathrm{eq}}(n_2,T_2)`$, we then observed in the simulations that $`S_{\mathrm{gas}}\left(f_1\right)`$ decrease while $`S_{\left(2\right)}(f,E)`$ increase during the evolution towards the equilibrium distributions.

Concluding remarks: We have confirmed via computer simulations the monotone increase of the Boltzmann entropy (log of phase space volume) for a dense fluid not in local equilibrium whose macrostate is specified by the empirical density $`f`$ and energy $`E`$. Similar results were obtained for a binary system of hard points in $`d=1`$. The simulations also show an apparent deterministic evolution of $`f_t`$ for such systems.

Acknowledgements: P.L.G acknowledges the support of MEC (Spain), the MCYT-FEDER project BFM2001-2841. The work of JLL was supported by NSF Grant DMR 01-279-26, and by AFOSR Grant AF 49620-01-1-0154. J.L.L. would like to thank O. Penrose for useful discussions.


text:
# Nonlocal potentials in independent-electron models

## I Introduction

The density functional theory (DFT) of Hohenberg, Kohn, and Sham has become the method of choice for computation of electronic structure and properties of large molecules and of solids. This theory justifies an independent-electron model that is considered to be exact in principle for N-electron ground states. It is widely assumed that such an exact model can be expressed in terms of local potential functions for exchange and correlation, although evidence to the contrary has existed for some time. It has only recently been shown that this locality hypothesis is inconsistent with exact theory for more than two electrons . The specific mathematical point is that density functional derivatives of the Hohenberg-Kohn universal functional cannot be equivalent to local potential functions, as assumed in established literature, for more than two electrons in a compact electronic system. These functional derivatives can be constructed explicitly, and are found to be Gâteaux derivatives, the generalization to functional analysis of analytic partial derivatives. For more than two electrons, they do not reduce as commonly assumed to simpler Fréchet derivatives, equivalent to multiplicative local potential functions, the generalization to functional analysis of analytic total derivatives.

This mathematical development explains a number of inconsistencies or paradoxes that result from assuming the locality hypothesis. Failure of locality for the exchange potential in a Hartree-Fock variational model of DFT is implied by results of optimized effective potential (OEP) calculations going back as far as Aashamar, Luke and Talman, confirmed by recent quantitative tests. The resulting problems or paradoxes include failure of DFT linear response theory to agree in the exchange-only limit with the time-dependent Hartree-Fock theory of Dirac. This failure can be traced to the locality hypothesis.

The most striking paradox affects the kinetic energy. In contrast to the locality hypothesis, Kohn and Sham use the kinetic energy operator of Schrödinger for what should be the density functional derivative of a ground-state kinetic energy functional. If the locality hypothesis were a rigorous consequence of variational theory, it should be valid for the kinetic energy functional, implying for noninteracting electrons an exact Thomas-Fermi theory (TFT) equivalent to DFT. However, the variational equations of these theories are incompatible for more than two electrons. The Kohn-Sham equations of DFT are correct for atomic shell structure, but the locality hypothesis in TFT is inconsistent with the exclusion principle. This is evident for the lowest $`1s2s^3S`$ state of an atom with two noninteracting electrons, perhaps the simplest unequivocal example of this inconsistency. Unless both $`1s`$ and $`2s`$ orbital densities are independently normalized, the lowest total energy is achieved for two electrons in the $`1s`$ state. Fermi-Dirac statistics require an independent Lagrange multiplier $`ϵ_i`$ for each orbital, while TFT provides only a single parameter $`\mu `$ for the total density. An equivalent inconsistency is implied for interacting electrons. This conclusion has been disputed and defended in recent publications.

Given a choice between an exact theory and one that can be expressed entirely in terms of local potential functions, there are pragmatic arguments supporting either alternative. The existing computational methodology of DFT requires use of local potential functions. In the local-density approximation (LDA), which is entirely justified by variational theory, first-principles calculations for condensed matter have had a revolutionary impact on understanding physical properties due to electronic structure. However, for magnetically ordered materials, there is a persistent problem due to the nonphysical electronic self-interaction, which may require use of nonlocal potentials in a fully predictive theory. For molecules, long-range potentials due to polarization response can become significant in determining molecular conformations, but apparently cannot be modeled using only approximate local correlation potentials. The optical potentials required in electron scattering theory are not directly described by energy-independent local potentials.

With such applications or extensions of existing theory in mind, the present paper will consider practical issues that arise in extending existing local-potential methodology to the nonlocal potential functions appropriate to an ultimately exact independent-electron model.

## II Variational theory

The variational theory of independent-electron models is most simply developed as an orbital functional theory (OFT). For an N-electron eigenstate such that $`(HE)\mathrm{\Psi }=0`$ and any rule $`\mathrm{\Psi }\mathrm{\Phi }`$ that determines a model or reference state $`\mathrm{\Phi }`$, unsymmetric normalization $`(\mathrm{\Phi }|\mathrm{\Psi })=(\mathrm{\Phi }|\mathrm{\Phi })=1`$ implies that $`E=(\mathrm{\Phi }|H|\mathrm{\Psi })=E_0+E_c`$, where $`E_0=(\mathrm{\Phi }|H|\mathrm{\Phi })`$ is an explicit orbital functional, and $`E_c=(\mathrm{\Phi }|H|\mathrm{\Psi }\mathrm{\Phi })`$ defines the correlation energy. Restricting the discussion to nondegenerate states, $`\mathrm{\Phi }`$ is a Slater determinant constructed from occupied orbital functions $`\{\varphi _i\}`$, with occupation numbers $`n_i=1`$. Spin indices and sums are assumed here but are suppressed in the notation. For orthonormal orbital functions, $`E_0=T+U+V`$, where

$`T={\displaystyle \underset{i}{}}n_i(i|\widehat{t}|i);V={\displaystyle \underset{i}{}}n_i(i|v|i),`$ (1)
for $`\widehat{t}=\frac{1}{2}^2`$. Two-electron functionals are defined for $`u=1/r_{12}`$ by $`U=E_h+E_x`$, where

$`E_h={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}n_in_j(ij|u|ij);E_x={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}n_in_j(ij|u|ji).`$ (2)
If $`𝒬=I\mathrm{\Phi }\mathrm{\Phi }^{}`$, $`E_c=(\mathrm{\Phi }|H|\mathrm{\Psi }\mathrm{\Phi })=(\mathrm{\Phi }|H|𝒬\mathrm{\Psi })`$. This implies an exact but implicit orbital functional

$`E_c=(\mathrm{\Phi }|H[𝒬(HE_0E_ci\eta )𝒬]^1H|\mathrm{\Phi }),`$ (3)
for $`\eta 0+`$.

Orbital Euler-Lagrange (OEL) equations follow immediately from standard variational theory, in terms of the orbital functional derivatives

$`{\displaystyle \frac{\delta T}{n_i\delta \varphi _i^{}}}=\widehat{t}\varphi _i;{\displaystyle \frac{\delta U}{n_i\delta \varphi _i^{}}}=\widehat{u}\varphi _i;`$ (4)
$`{\displaystyle \frac{\delta V}{n_i\delta \varphi _i^{}}}=v(𝐫)\varphi _i;{\displaystyle \frac{\delta E_c}{n_i\delta \varphi _i^{}}}=\widehat{v}_c\varphi _i,`$ (5)
using $`\widehat{u}=v_h(𝐫)+\widehat{v}_x`$, where $`v_h`$ is the classical Coulomb potential, and $`\widehat{v}_x`$ is the Fock exchange operator. For independent subshell normalization, Lagrange terms $`[(i|i)1]ϵ_i`$ are subtracted from the energy functional. Defining $`𝒢\varphi _i=\frac{\delta E}{n_i\delta \varphi _i^{}}`$, the variational equation is

$`{\displaystyle d^3𝐫\underset{i}{}n_i(\delta \varphi _i^{}\{𝒢ϵ_i\}\varphi _i+cc)}=0,`$ (6)
for unconstrained orbital variations in the usual Hilbert space. This implies the OEL equations

$`\{𝒢ϵ_i\}\varphi _i=0,i=1,\mathrm{},N.`$ (7)
Orbitals of different energy are orthogonal. Equivalently, the universal functional $`F=EV`$ defines

$`{\displaystyle \frac{\delta F}{n_i\delta \varphi _i^{}}}=\varphi _i=\{\widehat{t}+\widehat{u}+\widehat{v}_c\}\varphi _i.`$ (8)
and the OEL equations are

$`\varphi _i=\{ϵ_iv\}\varphi _i,i=1,\mathrm{},N.`$ (9)
Overlap products in $`(\mathrm{\Phi }|H|\mathrm{\Phi })`$ drop out of this argument, because they either produce a unit factor or are multiplied by a factor that vanishes. Explicitly, for noninteracting electrons, the same OEL equations are obtained if the simplified functional $`E_0=T+V`$ is replaced by $`(\mathrm{\Phi }|\widehat{T}+\widehat{V}|\mathrm{\Phi })`$. The functional differential

$`{\displaystyle \underset{i}{}}n_i`$ $`{\displaystyle d^3𝐫[\delta \varphi _i^{}(𝐫)\{\widehat{t}+v(𝐫)ϵ_i\}\varphi _i(𝐫)+cc]\mathrm{\Pi }_{ji}(j|j)}`$ (10)
$`+{\displaystyle \underset{i}{}}n_i`$ $`(i|\widehat{t}+v(𝐫)ϵ_i|i)\delta \mathrm{\Pi }_{ji}(j|j).`$ (11)
must vanish for unconstrained orbital variations about a stationary state. It cannot do so unless all Lagrange multipliers $`ϵ_i`$ are included, precluding an equivalent Thomas-Fermi theory. Variations of the overlap product are multiplied by a factor that vanishes.

## III The need for nonlocal potentials

Hohenberg-Kohn theory can be proven for a nondegenerate ground state of any version of OFT in which $`E_{xc}=E_x+E_c`$ is approximated by a parametrized functional of $`\rho `$. Model electronic density $`\rho =_in_i\rho _i`$ is a sum of orbital subshell densities $`\rho _i(𝐫)=\varphi _i^{}(𝐫)\varphi _i(𝐫)`$. Defining $`\widehat{v}_{xc}\varphi _i=\delta E_{xc}/n_i\delta \varphi _i^{}`$, the corresponding Kohn-Sham theory is determined by the OEL equations given above. The still controversial issue of whether or not the effective one-electron potentials in these equations can be expressed as multiplicative local potential functions depends on the relationship between the orbital functional derivatives of Eqs.(4) and the density functional derivatives needed for DFT.

Density functional derivatives must be consistent with the OEL equations. A crucial logical point is that because functional definitions must be extended to include variations that are not constrained by normalization, variational theory restricted to normalized ground states cannot determine OEL equations. A variational derivation of the Kohn-Sham equations must extend the universal ground-state functional of Hohenberg and Kohn to density variations driven by unrestricted variations in the orbital Hilbert space. This extended functional is denoted here by $`F_s`$. Orbital theory simplifies the mathematics, because such variations are necessary and sufficient to determine the Euler-Lagrange equations.

The rigorous theory of Englisch and Englisch proves the existence of a generic density functional derivative for a nondegenerate ground state. Although not considered in the original derivation, it is clear from the OEL Eqs.(9) that the undetermined constant in , Eq.(4.1), must have orbital indices in order to be consistent with the electronic subshell structure implied by Fermi-Dirac statistics. This is made explicit by Eqs.(8, 9) here, which imply that a functional derivative exists that depends on an orbital subshell index. Evaluated for a ground state, the functional differential

$`\delta F_s`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i{\displaystyle d^3𝐫\{\delta \varphi _i^{}(𝐫)\varphi _i(𝐫)+cc\}}`$ (12)
$`=`$ $`{\displaystyle \underset{i}{}}n_i{\displaystyle d^3𝐫\{ϵ_iv(𝐫)\}\delta \rho _i(𝐫)}`$ (13)
determines the functional derivative

$`{\displaystyle \frac{\delta F_s}{n_i\delta \rho _i}}=ϵ_iv(𝐫).`$ (14)
Because this depends on an orbital index, it defines a Gâteaux derivative. If a Fréchet derivative were to exist, equivalent to a multiplicative local potential function, it would imply the Thomas-Fermi equation,

$`{\displaystyle \frac{\delta F_s}{\delta \rho }}=\mu v(𝐫).`$ (15)
Because all $`\delta F_s/n_i\delta \rho _i`$ must equal $`\delta F_s/\delta \rho `$ if the latter exists, these equations are inconsistent unless all $`ϵ_i`$ are equal. Hence a Fréchet derivative cannot exist for more than one electron of each spin without violating the exclusion principle. This result is not in conflict with rigorous analysis, which establishes the existence of Gâteaux derivatives in general, but cannot distinguish between Gâteaux and Fréchet derivatives because normalization is constrained. For noninteracting electrons, this analysis applies to the Kohn-Sham kinetic energy functional $`T_s`$. If independent subshell normalization is enforced, using Gâteaux derivatives and the set of Lagrange multipliers $`\{ϵ_i\}`$, a generalized Thomas-Fermi theory can be expressed in terms of partial densities and is equivalent to the orbital functional theory of Eqs.(9).

The mathematical situation is that for more than one electron of either spin, the functional $`F_s`$ cannot be extended to derive Euler-Lagrange equations for an unstructured total density $`\rho `$, because a Fréchet functional derivative does not exist. However, because Gâteaux derivatives exist, an extended functional of orbital subshell densities does exist. Independent subshell variations imply Euler-Lagrange equations equivalent to the OEL equations, valid for electronic subshell structure and for the exclusion principle. For noninteracting electrons, $`F_s`$ reduces to the Kohn-Sham functional $`T_s`$, with the same mathematical properties. There is no implication that the density functional derivative of $`E_{xc}`$ must be a local potential function $`v_{xc}(𝐫)`$, unless this is mandated by the simplified functional dependence assumed in the local density approximation (LDA).

## IV Nonlocal potentials

The concept of a Gâteaux functional derivative establishes a direct relationship between orbital-indexed local potential functions and nonlocal potentials defined by orbital functional derivatives. Using exchange energy as an example, the orbital functional derivative is $`\delta E_x/n_i\delta \varphi _i^{}=\widehat{v}_x\varphi _i`$, and the functional differential can be expressed as

$`\delta E_x`$ $`=`$ $`{\displaystyle \underset{i}{}}n_i{\displaystyle d^3𝐫[\delta \varphi _i^{}\widehat{v}_x\varphi _i+cc]}`$ (16)
$`=`$ $`{\displaystyle \underset{i}{}}n_i{\displaystyle d^3𝐫\frac{\varphi _i^{}\widehat{v}_x\varphi _i}{\varphi _i^{}\varphi _i}\delta \rho _i}.`$ (17)
This determines a Gâteaux derivative

$`{\displaystyle \frac{\delta E_x}{n_i\delta \rho _i}}={\displaystyle \frac{\varphi _i^{}\widehat{v}_x\varphi _i}{\varphi _i^{}\varphi _i}}=v_{xi}(𝐫),`$ (18)
valid for arbitrary variations in the orbital Hilbert space. The technical problem of generalizing Kohn-Sham methodology to an exact theory reduces to methodology for such indexed local potentials. The occurrence of different potentials for different orbital indices can be compensated during self-consistent iterations by using off-diagonal Lagrange multipliers for orthogonalization. These parameters must vanish on convergence for closed-shell OEL equations. In the case of exchange in atoms, this is a variant of well-established methodology for iterative solution of the Hartree-Fock equations. Inclusion of an indexed correlation potential, derived from any parametrized functional $`E_c`$,

$`v_{ci}(𝐫)={\displaystyle \frac{\varphi _i^{}\widehat{v}_c\varphi _i}{\varphi _i^{}\varphi _i}}`$ (19)
is a straightforward generalization.

Hartree-Fock methodology based on expansion in atomic Gaussian basis orbitals is valid in principle, but becomes impractical for large molecules and solids. For solids, multiple scattering theory (MST) has been the method of choice for DFT calculations. It provides a common methodology for molecules, based on variational treatment of space-filling atomic cells. In MST, the principal calculation over a range of orbital energies is to solve the OEL equations within each local atomic cell, or within the enclosing sphere of each such cell. Independent representations of a wave function within each cell are matched across intercell boundaries using variational formalism. Generalization of the intracell calculations to nonlocal potentials is straightforward, simply including an indexed parametrized exchange-correlation potential while integrating Schrödinger or semirelativistic equations at specified energy. The practical technical problem is thus reduced to the treatment of the long-range tails of indexed local potentials defined within each cell. In MST, the formalism is greatly simplified by the fact that atomic basis functions are literally truncated at cell boundaries. This means that there are no intercell overlap contributions.

An approximation that may turn out to be both efficient and accurate is suggested by the MST formalism. This is a ”local nonlocal” (LNL) model. True indexed potentials are to be used only within atomic cells, but are replaced by their asymptotic forms in external cells. In general, these asymptotic forms arise from electrostatic multipole potentials, and can be represented as long-range local potentials. If this is valid, the indexed potentials in each cell are augmented by an external field potential that is summed over all other cells. As in the well-known Ewald expansion, this implies substantial cancellation of fields, including electrostatic screening, and may justify neglecting such terms during self-consistent iterations, including them as first-order perturbations after convergence.

An important step toward practical methodology for long-range correlation has been proposed, and tested in an electron scattering model. In a general context, the method is to compute, for each atomic cell, the first-order perturbed orbital functions due to external multipole fields. Intercell correlation energy then can be formulated in terms of these polarization pseudostates, and response potentials can be constructed that are correct at long range, but vanish smoothly at the cell origin. An example of dipolar response has been carried out in detail. The implied polarization potential is used to calculate polarization response effects in $`eHe`$ scattering cross sections.


text:
# Macroscopic evidence of quantum coherent oscillations of the total spin in the Mn-[3x3] molecular nanomagnet.

(October 2, 2003)

## Abstract

Molecular nanomagnets, besides promising to open new frontiers in technology, have attracted huge interest in the scientific community because they can exhibit the phenomenon known as quantum tunnelling of the magnetization, i.e. coherent fluctuations of the direction of the total spin vector. In this paper we study a different quantum phenomenon involving fluctuations of the magnitude of the total spin vector. These fluctuations are related to the mixing between states with different spin quantum number, and imply new macroscopic effects, which we theoretically investigated in the Mn-\[3x3\] grid.

Molecular nanomagnets (MNMs) Sessoli93 ; Gatteschi94 ; Legget95  are clusters containing a finite number of transition-metal ions whose magnetic moments (spins) are so strongly coupled that at low temperature each molecule behaves like a single-domain particle with fixed total spin. Being at the crossover between classical and quantum regimes, MNMs exhibit at the same time classical properties of macroscopic magnets such as magnetization hysteresis, and quantum phenomena like tunneling of the direction of the total spin through energy barriersFriedman96 ; Thomas96 ; Wernsdorfer99 . MNM systems are interesting also for potential technological applications, as envisaged for the implementation of quantum computing algorithms Leuenberger01 , or for dense and highly efficient memory devices Sessoli93 . Here we study a new macroscopic manifestation of a quantum phenomenon involving fluctuations not only of the direction, but also of the magnitude of the total spin, and we show that it is realized in a Mn-\[3$`\times `$3\] grid Thompson00 ; Wald02 . Recognizing the effects of these fluctuations is essential to achieve a satisfactory understanding of the role played by quantum mechanics in complex mesoscopic magnetic systems. Quantum magnetic phenomena were identified even in molecules of great biological interest such as ferritinferritin1 ; ferritin2 .

The advantage of studying quantum phenomena in MNMs is that the vanishingly small interaction between different molecules allows single-molecule phenomena to be observed at a macroscopic scale, because the crystal behaves like a collection of independent objects Gatteschi94 ; Wald97 , each described in general by spin Hamiltonians of the form Bencini90

$$H=\underset{i>j}{}J_{ij}𝐬_i𝐬_j+\underset{i}{}𝐬_i𝐃_i𝐬_i+$$

$$\underset{i>j}{}𝐬_i𝐃_{ij}𝐬_j+\mu _B\underset{i}{}g_i𝐁𝐬_i,$$
$`(1)`$
where $`𝐬_i`$ are spin operators of the $`i^{th}`$ ion in the molecule. The first term is the isotropic Heisenberg exchange interaction, the second and third terms describe the local crystal-field and the anisotropic intra-cluster spin-spin interactions. The last term is the Zeeman coupling with an external field $`𝐁`$ in which isotropic $`g`$-factors are assumed. While the Heisenberg term is rotationally invariant and therefore conserves the length $`|𝐒|`$ of the total spin $`𝐒=_i𝐬_i`$, the anisotropic terms break rotational invariance and do not conserve this observable. Nevertheless, since the Heisenberg contribution is usually largely dominant, $`|𝐒|`$ is nearly conserved, and the energy spectrum of $`H`$ consists of a series of level multiplets with an almost definite value of $`|𝐒|`$ (expressed in terms of the quantum number $`S`$ as $`\sqrt{S(S+1)}`$). Thus, quantum fluctuations of $`|𝐒|`$, which are associated with mixing of states with different value of the quantum number $`S`$ (”$`S`$-mixing” Liviotti02 ), either are zero or are expected to produce negligible effects on the macroscopic behavior, and are therefore neglected in virtually all studies. A major theoretical goal would be to identify a clear macroscopic signature of such fluctuations, and a model system displaying such effect.

In this paper we show that favorable conditions are met in the Mn(II)-\[3$`\times `$3\] grid-like cluster \[Mn9(2POAP-2H)6\](ClO4)6$``$3.57MeCN$``$H2O (hereafter Mn-\[3$`\times `$3\], see Fig. 1), a system that has been recently characterized by magnetization and torque measurements Wald02 ; torque . Our theoretical calculations suggest grid-shaped molecules as good candidates to study fluctuations of $`|𝐒|`$, since the lowest level is expected to display significant $`S`$-mixing with the first excited multiplet. In fact, unlike in ideal ring-shaped molecules Waldepl , in the \[3$`\times `$3\] grid the two lowest manifolds of the exchange part of the spin Hamiltonian belong to the same irreducible representation of the molecular symmetry group. Therefore, the application of a suitably oriented magnetic field induces a series of anticrossings (ACs) between the ground state and levels originating from higher excited manifolds (see red lines in Fig. 2). At the AC fields, $`S`$-mixing in the ground state is maximum, and quantum fluctuations of $`|𝐒|`$ are greatly enhanced. TorqueWald02  and neutronGuidiUn  experiments on Mn-\[3$`\times `$3\] show that the zero-field gap between the two lowest $`S`$-multiplets is small enough for the ACs to occur at fields within experimental reach.

Mn-\[3$`\times `$3\] crystallizes in the space group $`C_2/c`$, and the cation \[Mn9(2POAP-2H)6\]6+ exhibits a slightly distorted $`S_4`$ molecular symmetry with the $`C_2`$ axis perpendicular to the plane of the cluster Thompson00 . The average distance between the Mn(II) ions is 3.93 $`\AA `$, the smallest distance between clusters is larger than 8 $`\AA `$. For a cluster composed of nine interacting Mn(II) spins with $`𝐬_i=5/2`$ the dimension of the Hilbert space is 10077696. The difficulties related with this huge dimension have been overcome by exploiting both the irreducible tensor operator technique and the Lanczos algorithm for the exact diagonalisation. The two-step procedure already developed Carretta03  has allowed the inclusion of $`S`$-mixing effects in the calculation. The exchange integrals and the single-site and spin-spin anisotropy tensors have been determined by neutron spectroscopy GuidiUn . All exchange integrals between nearest-neighbors are found to be nearly equal to $`0.47`$ meV, apart from $`J_{18}`$, $`J_{78}`$, $`J_{34}`$ and $`J_{45}`$ (see Fig. 1), whose values are $`0.33`$ meV. Next-nearest-neighbor exchange interactions can be neglected. Concerning the local crystal-field, the second term in Eq. 1 can be approximately rewritten as

$$\underset{i}{}𝐬_i𝐃_i𝐬_i=D\underset{i}{}[s_{iz}^2\frac{1}{3}s_i(s_i+1)],$$
$`(2)`$
with $`D=6.1`$ $`\mu `$eV. The intracluster dipole-dipole interaction $`𝐃_{ij}`$ has been evaluated within the point-dipole approximationBencini90 . At last, $`g_i=2`$ is assumed, as appropriate for Mn(II) ions.

With these experimentally determined parameters, the energy spectrum for $`𝐁=0`$ consists of many level multiplets with an almost definite value of $`S`$. The ground multiplet has $`S=5/2`$, and the four lowest-lying excited multiplets have (in order of increasing energy) $`S=7/2`$, $`S=3/2`$, $`S=3/2`$, $`S=9/2`$. These multiplets are separated by the isotropic exchange and split by the anisotropic interactions. The lowest levels are shown in Fig. 2. The application of a magnetic field $`𝐁`$ in a direction outside the grid plane and different from that of the $`C_2`$ axis produces several ACs involving levels belonging to different multiplets (see Fig. 2). As the AC fields $`𝐁_𝐜`$ are approached the multiplet mixing is enhanced. For $`𝐁=𝐁_𝐜`$, the spins in each cluster oscillate coherently between states with different values of $`S`$, which therefore is no longer a good quantum number. This can be inferred for example for the ground state ACs (indicated by the arrows in Fig. 2) from the inset in Fig. 2, which shows the field-dependence at $`T=0`$ of the quantity $`S_{eff}`$, defined through the relation $`S_{eff}(S_{eff}+1)=𝐒^2`$. When $`𝐁`$ is close to $`𝐁_𝐜`$ the value of $`S_{eff}`$ is intermediate between half integers, e.g. 5/2 and 7/2 at the first level AC, thus confirming that the ground-state wavefunction is a superposition of different total-spin states.

This opens a scenario in which not only the direction of the total spin fluctuates in time, as in the case of the well-studied quantum tunnelling of the magnetization Friedman96 ; Thomas96 ; Wernsdorfer99 , but even its length fluctuates (quantum dynamics of the total spin, or briefly QDTS). Moreover, by properly tuning the direction of the applied magnetic field, the AC splitting can be made as large as several Kelvins, thus overcoming the problem of decoherence due to hyperfine fields and to cluster-cluster interactions. Most importantly, in the Mn-\[3$`\times `$3\] grid several of these ACs involve the two lowest energy levels, thus opening the possibility to detect the QDTS by means of $`\mathrm{𝑚𝑎𝑐𝑟𝑜𝑠𝑐𝑜𝑝𝑖𝑐}`$ low-temperature magnetic bulk techniques. In particular, our calculations show that at low temperature, in correspondence to each AC between the two lowest lying states, a sharp peak should be detected in the torque signal (as function of the field intensity $`B`$) when the direction of the field $`\theta `$ is close to the grid plane ($`\theta =0^o`$) or almost perpendicular to it ($`\theta =90^o`$) (see Figs. 3 and 4). The precipitous drop of the torque signal as a function of $`\theta `$ near $`\theta =0^o`$ and $`\theta =90^o`$ is due to a change of sign imposed by symmetry (as $`𝐌\times 𝐁=0`$ when $`𝐁`$ points along a principal direction of the susceptibility tensor of the system).

By comparing Figs. 3 and 4 with the corresponding spectrum reported in Fig. 2, it is evident that a torque peak appears in correspondence to each AC involving the ground state. While the first peak at low field corresponds to an intra-multiplet AC, the peaks at about $`6.6`$, $`8.9`$ and $`11.2`$ Tesla represent a direct consequence of the coherent superposition of different total spin quantum states. The link between the peaks in the torque signal and the level ACs is demonstrated in Fig. 4, where the torque curves calculated with the mixing artificially forced to zero are shown by colored dashed lines. When $`S`$-mixing is neglected, every inter-multiplet AC becomes a crossing, the fluctuations of the total spin vector are suppressed and the peaks in the torque curve disappear. Thus, without $`S`$-mixing the torque curve should exhibit the step-like field dependence usually observed in MNMs Waldtorque . The comparison between solid and dashed curves in Fig. 4 shows clearly that the predicted peaks in the torque would be a direct macroscopic manifestation of the QDTS at the level ACs. Also, these spin fluctuations and the associated peaks are of purely quantum origin, as these are absent in the classical version of the Hamiltonian Eq. (1) (see the gray dashed curve in Fig. 4).

The fact that these torque peaks directly reflect the large enhancement of quantum fluctuations of $`|𝐒|`$ at the ACs can be understood by a simple physical picture. In the following we define the $`z^{}`$ axis parallel to $`𝐁`$. The torque is proportional to the magnetic response perpendicular to the field direction, i.e. to $`S_x^{}`$. Fluctuations of $`|𝐒|`$ are accompanied by fluctuations of the magnetization ($`S_z^{}`$)collinear , and these latter are deeply connected with $`S_x^{}`$. Accordingly, $`S_x^{}`$ tracks the increase and decrease of these fluctuations while sweeping over the AC, leading to a peak in the torque. Indeed, near the ACs, the ground-state wavefunction can be written to a good degree of approximation as

$$|G=\frac{1}{\sqrt{2}}(a(B)|\mathrm{\Gamma }_1,M+b(B)|\mathrm{\Gamma }_2,M+1),$$
$`(3)`$
where $`a(B)^2+b(B)^2=2`$ by normalization. For $`BB_c`$, $`b(B)0`$, $`a(B)\sqrt{2}`$, and the ground state reduces approximately to an eigenstate of $`S_z^{}`$, $`|\mathrm{\Gamma }_1,M`$, where $`M`$ is the corresponding eigenvalue and $`\mathrm{\Gamma }_1`$ represents the set of additional labels necessary to identify the state. For $`BB_c`$, $`a(B)0`$, $`b(B)\sqrt{2}`$, and the ground state reduces approximately to $`|\mathrm{\Gamma }_2,M+1`$.

At $`T=0`$ the torque is given by

$$\tau BS_x^{}=Ba(B)b(B)\mathrm{\Gamma }_1,M|S_x^{}|\mathrm{\Gamma }_2,M+1.$$
$`(4)`$
Since $`\mathrm{\Gamma }_1,M|S_x^{}|\mathrm{\Gamma }_2,M+1`$ is almost independent on $`B`$, the field dependence of $`\tau `$ comes entirely from the factor $`Ba(B)b(B)=Bb(B)\sqrt{2b(B)^2}`$.

Quantum fluctuations of the magnetization are usually characterized by the quantity

$$(\mathrm{\Delta }S_z^{})^2=S_z^{}^2S_z^{}^2=\frac{b(B)^2}{2}\frac{b(B)^4}{4}.$$
$`(5)`$
Therefore,

$$\tau 2B\mathrm{\Delta }S_z^{}.$$
$`(6)`$
In case of anticrossing, $`\mathrm{\Delta }S_z^{}`$ is maximum at $`B_c`$, where $`b(B_c)=a(B_c)=1`$, and the torque peaks. In case of crossing $`\mathrm{\Delta }S_z^{}`$ is always zero because either $`b(B)=0`$ or $`a(B)=0`$, and the torque does not peak. As a further confirmation of this picture we show in Fig. 5 numerical calculations of the fluctuations of $`S_z^{}`$ (”parallel”) and of $`S_x^{}`$ (” perpendicular”), with and without $`Smixing`$.

The torque is the macroscopic magnetic bulk technique displaying the clearest signature of $`S`$-mixing. When the temperature $`T`$ increases these effects wash out, and indeed published torque measurements on Mn-\[3$`\times `$3\] Wald02  do not display peaks because data were collected at too high a temperature.

We have also calculated the low-$`T`$ heat capacity as a function of the applied field (see Fig. 6). This quantity provides a way to assess the value of the AC gaps as a function of the field direction. For $`\theta 0`$ anticrossings open up and the heat capacity does not vanish at $`B_c`$.

In this paper we have identified a new clear-cut macroscopic manifestation of a quantum phenomenon characterizing the microscopic dynamics of magnetic clusters. We have shown that torque measurements provide direct evidence of quantum fluctuations of the total spin length, and we have studied a model system displaying such effect. The comprehensive understanding of the spin dynamics in mesoscopic systems like molecular nanomagnets is essential to understand the crossover between quantum and classical mechanics and to make rational design of these compounds possible, especially in view of the promising applications.

Fig. 1: A schematic representation of the molecular structure of the Mn-\[3$`\times `$3\] grid showing the network of manganese (blue) and oxygen (green) bonds. Other atoms are omitted for clarity. The red arrows indicate the spin directions in the classical ground state.

Fig. 2: Calculated energy levels of the Hamiltonian Eq.(1) with the parameters given in the text, and with the direction of the applied field forming an angle $`\theta =2.8^o`$ with the grid plane. Energies are plotted as functions of the applied field intensity $`B`$. The ground state energy is set equal to zero. Arrows indicate the anticrossings produced by $`S`$-mixing between the two lowest levels (highlighted in red). The inset shows $`S_{eff}`$ as functions of $`B`$ for the same angle at $`T=0`$, where $`S_{eff}`$ is defined through the relation $`S_{eff}(S_{eff}+1)=𝐒^2`$.

Fig. 3: Two-dimensional plot of the calculated torque as a function of the applied field direction $`\theta `$ and intensity $`B`$ at $`T=0.05`$ K with parameters of Hamiltonian Eq.(1) as given in the text.

Fig. 4: Plots of the calculated torque at $`T=0.4`$ K as a function of the applied magnetic field $`B`$ for two selected directions $`\theta `$ lying close to the grid plane. Colored dashed lines represent the calculated torque at the same angles with $`S`$-mixing eliminated by hand, i.e. with quantum fluctuations of $`|𝐒|`$ removed. The dashed gray line represents the torque calculated using the classical version of the Hamiltonian Eq. (1) at $`T=0`$ for $`\theta =2.8^o`$. The classical value of the torque has been rescaled (by about 3) to fit the figure. The inset shows the calculated torque at $`T=0.4`$ K as a function of the applied field intensity $`B`$ for several directions $`\theta `$.

Fig. 5: Numerical calculations of the $`T=0`$ fluctuations of $`S_z^{}`$ (”parallel”) and of $`S_x^{}`$ (” perpendicular”), as a function of the applied magnetic field (forming an angle $`\theta =2.8^o`$ with the grid plane), with and without $`Smixing`$. Parameters are given in the text.

Fig. 6: Numerical calculations of the heat capacity as a function of the applied magnetic field for several directions $`\theta `$, at $`T=0.4`$ K.