text:
# Relaxation, the Boltzmann-Jeans Conjecture and Chaos

## Abstract

Slow (logarithmic) relaxation from a highly excited state is studied in a Hamiltonian system with many degrees of freedom. The relaxation time is shown to increase as the exponential of the square root of the energy of excitation, in agreement with the Boltzmann-Jeans conjecture, while it is found to be inversely proportional to residual Kolmogorov-Sinai entropy, introduced in this Letter. The increase of the thermodynamic entropy through this relaxation process is found to be proportional to this quantity.

Study of the relaxation from highly excited states is important not only in physico-chemical systems but also in biological systems. It has been reported that excitations of some protein molecules are maintained over anomalously long time spans . Such behavior is relevant to enzymatic reactions, and in particular to their biological functions, because these molecules must maintain their excited state (brought about, say, by ATP) to be able to proceed from one process to the next.

Typically, a Hamiltonian system with a sufficiently large number of degrees of freedom relaxes to equilibrium. This relaxation process, however, is not necessarily fast. For instance, some Hamiltonian systems can exhibit logarithmically slow relaxation from excited states, called Boltzmann-Jeans conjecture (BJC), which was first noted by Boltzmann, explored by Jeans, Landau and Teller, and then has been studied in terms of nonlinear dynamics . In the BJ conjecture, the relaxation to equilibrium is required, but existence of chaos is not explicitly assumed. On the other hand, irreversible relaxation is often studied in relationship with chaos. In the present paper we study slow relaxation of a type in agreement with the BJC, by using a Hamiltonian dynamical system with a large number of degrees of freedom, and explore its possible relationship with chaos.

We consider dynamics with the Hamiltonian

$$H=K+V=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2}+\underset{i,j=1}{\overset{N}{}}V(\theta _i,\theta _j)$$
(1)
describing a set of pendula, where $`p_i`$ is the momentum of the $`i`$-th pendulum. The potential $`V`$ is chosen so that each pair of pendula interacts only through their phase difference :

$$V(\theta _i,\theta _j)=\frac{1}{2(2\pi )^2N}\{1\mathrm{cos}(2\pi (\theta _i\theta _j))\}.$$
(2)
Hence, the evolution equations for the momentum $`p_i`$ and the phase $`\theta _i`$ are given by

$$\dot{p}_j=\frac{1}{2\pi N}\underset{i=1}{\overset{N}{}}\mathrm{sin}(2\pi (\theta _i\theta _j)),\dot{\theta }_j=p_j.$$
(3)
This form of the interaction is chosen so that there is an attractive force tending to align the phases of all the pendula. The thermodynamic properties of this model in the $`N\mathrm{}`$, as well as the finite-size effect, are investigated in Ref.. Here we mainly discuss the case with $`N=10`$, which is sufficiently large to exhibit the thermodynamic properties, and the results to be discussed are valid for a larger system. The temperature of the system can be defined as $`T2K/N`$, which is a monotonically increasing function of the total energy $`E`$.

As studied in Refs., each pendulum in this system exhibits small-amplitude vibration when the total energy is small, while, as the energy is increased, some pendula begin to display rotational behavior over many cycles. The relaxation of the rapidly rotating pendula into a vibrating assembly is rather slow, since their average interaction with the assembly almost cancels out over the slow time scale of the assembly, due to the rapid rotation. Here we concentrate on such slow relaxation of the rotational mode to the vibrational modes.

The slow relaxation from a highly excited state is investigated systematically in the situation that the excited state is prepared by applying an instantaneous kick to a certain pendulum in the system. An example of the relaxation process following this kick is depicted in Fig.1, where it is seen that the kicked pendulum continues to rotate in isolation, maintaining a large energy. The relaxation time $`t_R`$ is defined as the interval of time from the kick to the point of which kinetic energy of the kicked pendulum first becomes smaller than $`K/N`$.

The relaxation process depends on the energy $`E_0`$ imported by the kick and also on the state of the other $`N1`$ pendula. Hereafter we call this assembly the ‘bulk’. The bulk state is parameterized by the total energy $`E`$ (or the temperature $`T`$) of the system before the kick .

The mean relaxation time $`t_R`$ over an ensemble of initial conditions was computed by fixing $`E`$ and varying $`E_0`$. As shown in Fig.2, the relaxation time increases exponentially with the kicked energy $`E_0`$  as

$$t_R\mathrm{exp}(\alpha \sqrt{E_0}),$$
(4)
for sufficiently large values of $`E_0`$. Here $``$ represents an ensemble average and $`\alpha `$ is a constant. When the value of $`E_0`$ is fixed, $`t_R`$ decreases as a function of $`E`$, as is seen in the figure.

The above exponential form of the relaxation agrees with that given by the BJC, which describes the relaxation time for a system consisting of two parts with very different frequencies as an exponential of the ratio of the two frequencies. Assuming the separation of time scales of the kicked pendulum and the bulk pendula, this exponential growth can be obtained through the Fourier analysis of the interaction term, as shown by Landau and Teller  (see also ).

In this argument, chaotic dynamics are not explicitly included, while chaos is often relevant to the irreversible relaxation to equilibrium. Actually, we have observed the exponential law (4) in numerical studies clearly, as long as the Hamiltonian dynamics for the given energy $`E`$ uniformly exhibits chaotic behavior without remnants of KAM tori. In this regime, $`t_R`$ decreases with the energy $`E`$, as shown in the inset of Fig.2. On the other hand, in the regime of the small values of $`E`$, the relaxation process highly depends on the initial condition and it is difficult to define average relaxation time. The hindrance of relaxation is associated with the weakness of chaos, since the relaxation of the rotating pendulum stops when the dynamics of the bulk is trapped to the vicinity of remnants of tori. Now, we study a possible relationship between the relaxation time and the chaotic dynamics of the bulk.

First let us consider the nature of the dynamics in the limit $`E_0\mathrm{}`$. In this case, the kicked pendulum is completely free from the bulk of vibrating pendula, whose collective motion is simply that determined by the Hamiltonian of its $`(N1)`$ pendula. At a finite value of $`E_0`$, the dynamics of the bulk deviate from those in the $`E_0\mathrm{}`$ limit, due to the interaction with the kicked pendulum. In order to quantify such deviation, we choose bulk states with energies distributed around a given $`E`$, and compute the Lyapunov spectrum $`\lambda _i(E_0)`$ for various values of $`E_0`$. Note that, in the limit $`E_0\mathrm{}`$, $`\lambda _{N2}=0`$ corresponding to the free rotation of the kicked pendulum, while $`\lambda _{N1}=\lambda _N=0`$ always holds due to the conservation of the total momentum and energy.

Since the kinetic energy $`K_r`$ of the kicked pendulum fluctuates and relaxes, $`\lambda _i(E_0)`$ is computed in the case that $`K_r`$ satisfies $`|E_0K_r|<\delta `$, with a small constant $`\delta `$. Since the relaxation of $`K_r`$ is very slow, this condition is satisfied over a time sufficiently long to compute the local exponents. In order to insure the convergence of the Lyapunov exponents, we take an ensemble average of such local exponents . (The existence of such time scale for the convergence is guaranteed by the slow relaxation, and indeed is confirmed numerically.)

The values of the computed Lyapunov exponents increase as $`E_0`$ decreases. The increased part of the Lyapunov exponents from those at $`E_0\mathrm{}`$ limit characterizes the amplification of the chaotic instability due to the interaction between the kicked pendulum and the bulk. With this in mind, we define the residual Lyapunov spectra by

$$\mathrm{\Delta }\lambda _i(E_0)\lambda _i(E_0)\underset{E_0\mathrm{}}{lim}\lambda _i(E_0).$$
(5)
¿From detailed numerical experiments, we find that the two scaling relationships;

$`\lambda _{N2}(E_0)=\mathrm{\Delta }\lambda _{N2}(E_0)\mathrm{exp}(\beta _1\sqrt{E_0}),`$ (6)
$`\mathrm{\Delta }h(E_0){\displaystyle \underset{i=1}{\overset{N3}{}}}\mathrm{\Delta }\lambda _i\mathrm{exp}(\beta _2\sqrt{E_0})`$ (7)
hold for sufficiently large values of $`E_0`$. Here $`\beta _1`$ is smaller than $`\beta _2`$, as shown in Fig.3(a).

Since $`\lambda _{N2}=0`$ in the limit $`E_0\mathrm{}`$, the scaling of $`\lambda _{N2}`$ represents the decrease in the degree of chaos of the kicked pendulum. The sum of the remaining positive Lyapunov exponents, then, is expected to correspond to the Kolmogorov-Sinai (KS) entropy of the bulk, and the increase in the degree of chaos of the bulk due to the interaction with the kicked pendulum is given by $`\mathrm{\Delta }h`$. (Note that the increase of Lyapunov exponents is not due to the increase of the average velocity in the bulk. Indeed, we have numerically confirmed that the temperature of the bulk for finite $`E_0`$ is not shifted from that in the limit $`E_0\mathrm{}`$.)

Since $`\mathrm{\Delta }h`$ is associated with the dynamics of the bulk and how it is affected by the energy absorbed from the kicked pendulum, we expect that it is related to the relaxation time $`t_R`$ . Indeed $`\alpha `$ in Eq.(4) and $`\beta _2`$ in Eq.(7) appears to be equal from our numerical results. The relationship between $`t_R`$ and $`\mathrm{\Delta }h`$ is plotted in Fig.3(b), which supports the relationship

$$\mathrm{\Delta }ht_R^1.$$
(8)
Thus, we find that the degree of chaos is closely related to the relaxation time scale.

The relationship between the relaxation and the residual KS entropy $`\mathrm{\Delta }h`$ can be roughly explained as follows: In a Hamiltonian system, for every orbit, there exists a time-reversed orbit, obtained by changing the sign of all $`p_i`$. However, because of the interaction with the kicked element, to obtain the time-reversed orbit of the bulk, a slight modification of each pendulum, in addition to the simple reversed of each of their momenta, $`p_ip_i`$, is needed. Since the KS entropy $`h`$ is a measure of the time rate of loss of the information concerning the initial conditions as an orbit evolves, the measure of orbits that can be considered the reversed orbit of some given orbit over a time $`t`$ decreases with $`\mathrm{exp}(ht)`$. Therefore, the time scale for the irreversible relaxation is proportional to the inverse of the KS entropy difference between the original and the reversed orbits. Hence, the relationship Eq.(8) between the time scale for the relaxation and the residual KS entropy is reasonable.

Next, we study dynamics on the long time scale of the order of $`t_R`$. To compare the relaxation from various kicked energies $`E_0`$, it is convenient to use the scaled time $`\tau =t/t_R`$, where $`t_R`$ is a function of $`E_0`$. We note that the time scale of relaxation decreases as the energy of the kicked pendulum, $`K_r`$, decreases. Then the above definition of $`\tau `$ implies $`\tau =t/t_R(\stackrel{~}{K}_r)`$ along the relaxation, where $`\stackrel{~}{K}_r`$ is the coarse-grained energy of the kicked pendulum over a long time scale and $`t_R\mathrm{exp}(\alpha \sqrt{\stackrel{~}{K}_r})`$.

Now, we derive an equation describing the relaxation of $`\stackrel{~}{K}_r`$ on the time scale of $`\tau `$. We have computed the difference $`\delta K_r`$ over the time steps $`\delta \tau `$ for many samples of orbits starting from a given $`K_r`$. As shown in Fig.4, the data are satisfactory fit by $`\frac{\delta K_r}{\delta \tau }=C,`$ with a constant $`C`$ ($`0.035`$) for various values of $`K_r`$ and $`E`$. Then the coarse-grained equation for the relaxation is obtained as $`dK_r/d\tau =C.`$ Noting that the interaction energy between the bulk and the kicked pendulum is tiny and recalling the conservation of energy, we get $`dE/d\tau =C.`$ Thus, with regard to the relaxation process, the difference in the details of the dynamical properties of the system under different conditions are eliminated by considering the scaled time $`\tau `$, while the dynamics and properties of the bulk seem to be strongly dependent on the value of $`E`$. ¿From the above equation, for the original time variable $`t`$, $`K_r`$ decays in the logarithmic time scale.

Now it is possible to consider the entropy increase of the bulk occurring during the relaxation. Since we have assumed that the slow change of the bulk part can be described by the slow change of thermodynamic quantities, the increase of the entropy of the bulk is estimated as

$$\frac{dS}{d\tau }=\frac{dS}{dE}\frac{dE}{d\tau }=\frac{C}{T},$$
(9)
where $`T`$ is the temperature of the bulk. This relation is obtained by noting that the process here is sufficiently slow to define these thermodynamic quantities on a coarse-grained time scale.

Since $`\tau =t/t_R\mathrm{\Delta }ht`$, by Eq.(8), the relation

$$\frac{dS}{dt}\mathrm{\Delta }h$$
(10)
for the relaxation process is obtained from Eq.(9). The proportionality in Eq. (10) relates two kinds of directed motion towards equilibrium: $`\mathrm{\Delta }h`$ corresponds to the amplification of chaos in the equi-energy surface of the bulk, whereas $`dS/dt`$ represents the growth of this equi-energy surface towards equilibrium. It should also be noted that this relationship holds not near the equilibrium state of the total system but for the nonlinear relaxation process from a highly excited state.

The relationship between chaotic dynamics and irreversibility has been extensively studied . A relationship between the thermodynamic entropy and the irreversible information loss was proposed by Sasa and Komatsu in the case that an external operation causes a transition from one equilibrium state to another. We believe that our relation between the residual KS entropy and the thermodynamic entropy is related to it, although ours applies to the spontaneous relaxation.

To sum up, we have obtained the following relationship between the relaxation time, $`t_R`$, and the energy supplied by an external kick, $`E_0`$: $`t_R\mathrm{exp}(\alpha \sqrt{E_0})`$. Then, this relaxation time $`t_R`$ is found to be proportional to the inverse of the residual KS entropy, that is, the difference between the KS entropy of the bulk interacting with the kicked element and that of the isolated bulk. By rescaling time by $`t_R`$, the excited energy is found to relax at a rate that is independent of the energy of the kicked element and the bulk. Finally, we find evidence that the rate of entropy generation is proportional to the residual KS entropy.

Although these three relationships were found through numerical study of the Hamiltonian (1), we expect that they hold generally for relaxation processes in Hamiltonian systems under the condition that in the limit of high excitation, the interaction between the excited elements and the unexcited bulk vanishes. Indeed, we have numerically found that Eqs.(4) and (8) are valid also for a coupled pendulum model with the nearest-neighbor coupling on a square lattice. It is important to examine the universality of the three relationships we have found, and also to study their relevance in regard to intra-molecular relaxation processes, including application to the energy transduction carried out by biological molecular motors.

The authors are grateful to S. Sasa and T. S. Komatsu for stimulating discussions. This research was supported by Grants-in-Aids for Scientific Research from JSPS and the REIMEI Research Resources of JAERI.


format: text

text:
# Kinetic Theory of Traffic Flows

## I Introduction

Traffic flows are strongly interacting many-body systems. Therefore, theoretical techniques such as kinetic theory and hydrodynamics are useful in describing the rich phenomenology of traffic flows which includes shock waves, phase transitions, clustering, metastability, hysteresis, etc. Traffic is typically modeled within macroscopic descriptions such as hydrodynamics and kinetic theory  or microscopic approaches, e.g., cellular automata  and car-following models . The large body of recent work on the physics of traffic flows is surveyed in Refs. .

In this review, we describe how quenched randomness in the car velocities leads to formation of clusters in one lane roadways. We assume ballistic motion with infinite memory, namely that each car has a preferred “intrinsic” velocity by which it drives in the absence of other cars . While the emerging behavior is quite similar to that found in stochastic particle hopping processes with quenched disorder , these ballistic motion models often admit deeper analytical treatment. Our starting point is an idealized no passing flow, where an exact analytical solution is possible, and an exact kinetic theory can be constructed . We then treat more realistic generalizations where passing is allowed using approximate kinetic theories . Our goal is to provide a concise summary where key features are emphasized and outstanding open issues are highlighted.

## II No passing zones

Our basic traffic model mimics cluster formation (often also called platoon formation) in one-lane roadways where passing is forbidden . In this model, each car moves ballistically at its initial velocity until it overtakes the preceding car or cluster. After this encounter, the incident car assumes the velocity of the cluster which it has just joined. Cars are taken to be size-less, and collisions to be instantaneous. We primarily consider spatially homogeneous situations where the positions and the velocities of the cars are initially uncorrelated. Specifically, cars are distributed randomly in space with a concentration $`c_0`$, and their velocities are drawn from the initial velocity distribution $`P_0(v)`$. Remarkably, analytic expressions can be obtained for the velocity distribution and the joint size-velocity distribution of clusters. Furthermore, it is also possible to describe analytically spatial inhomogeneities and even input of cars into the system.

### A The velocity distribution

We first consider the cluster velocity distribution. In this description, only the lead car is relevant and trailing cars in a cluster can be ignored. Let $`P(v,t)`$ be the distribution of clusters with velocity $`v`$ at time $`t`$. Initially, all cars are lead cars and the cluster (or lead-car) velocity distribution equals $`P(v,t)=P_0(v)S(v,t)`$, a product of the initial velocity distribution and the survival probability $`S(v,t)`$. The survival probability is the probability that a car with initial velocity $`v`$ avoids “collisions” with slower cars up to time $`t`$, and hence, is still moving at the same velocity. Consider a car of initial velocity $`v`$. To ensure that it would not overtake slower cars of velocity $`v^{}<v`$, an interval of size $`(vv^{})t`$ ahead of it must not contain $`v^{}`$-cars initially. For the velocity distribution, $`P_0(v)`$, and a Poissonian initial spatial distribution, the probability for such an event is $`\mathrm{exp}\left[t(vv^{})P_0(v^{})\right]`$. For a car to survive to time $`t`$, this exclusion probability should be taken into account for every $`v^{}<v`$, and taking the product over all possible slower cars yields the survival probability, $`S(v,t)=\mathrm{exp}[t_0^v𝑑v^{}(vv^{})P_0(v^{})]`$. Consequently, the cluster velocity distribution is found for arbitrary initial conditions

$$P(v,t)=P_0(v)\mathrm{exp}\left[t_0^v𝑑v^{}(vv^{})P_0(v^{})\right].$$
(1)
The process is deterministic (the initial condition is the only source of randomness) and given the initial positions and velocities of the cars, the state of the system at any later time follows. This is reflected nicely in Eq. (1).

The exact solution (1) satisfies the following Boltzmann equation,

$$\frac{P(v,t)}{t}=P(v,t)_0^v𝑑v^{}(vv^{})P(v^{},0).$$
(2)
Interestingly, this rate equation is linear in the velocity distribution $`P(v,t)`$. The collision rate is proportional to the relative velocity, $`vv^{}`$, and the initial velocity distribution of the slower cars, a signature of the infinite memory in the system. This is a unique case where the hierarchy of evolution equations corresponding to the velocity distributions terminates at the first order.

Whereas the steady state of one-lane traffic with no passing is trivial, viz. all cars will eventually join a cluster led by the slowest car in the system, the time dependent behavior is interesting. We concentrate on the long time behavior, which is largely independent of the initial distribution of fast cars, as follows from Eq. (1). For discrete velocity distributions, the time dependence is exponential, and we focus on continuous distributions. In this case, we find directly from Eq. (1) that both the cluster concentration, $`c(t)=𝑑vP(v,t)`$, and the average velocity, $`v(t)=c^1𝑑vvP(v,t)`$, decay algebraically with time

$`c(t)`$ $``$ $`t^\alpha \alpha ={\displaystyle \frac{\mu +1}{\mu +2}},`$ (3)
$`v(t)`$ $``$ $`t^\beta \beta ={\displaystyle \frac{1}{\mu +2}}.`$ (4)
The scaling exponents $`\alpha `$ and $`\beta `$ depend only on the small-$`v`$ extremal statistics  of the initial velocity distribution,

$$P_0(v)av^\mu v0,$$
(5)
via the cutoff exponent $`\mu >1`$. The two scaling exponents are related by $`\alpha +\beta =1`$ as dictated by an elementary mean free path argument: $`cvt1`$.

Since the number of cars is conserved, the average cluster size is inversely proportional to the concentration, $`m1/c`$, and the size growth law is $`mt^\alpha `$. In the limit $`\mu \mathrm{}`$, the size grows linearly with time. In contrast, when $`\mu 1`$, the size remains roughly constant, since the velocity distribution becomes effectively unimodal and collisions become exceedingly rare. This qualitative dependence on the form of the initial velocity distribution is reminiscent of ballistic annihilation processes, where ballistically moving particles annihilate upon collision . The above clustering process can be viewed as a ballistic aggregation process that possesses a single mass conservation law. The sensitive dependence on the initial conditions is in contrast with momentum conserving ballistic agglomeration processes (that mimics large scale formation of matter in the universe) where a universal scaling asymptotic behavior emerges .

The average velocity is the only relevant velocity scale in the problem and asymptotically the velocity distribution follows the scaling form

$$P(v,t)t^{\beta \alpha }\mathrm{\Phi }(vt^\beta ).$$
(6)
From Eqs. (1) and (5), the scaled distribution is $`\mathrm{\Phi }(z)=az^\mu \mathrm{exp}(bz^{\mu +2})`$, with $`b=a/[(\mu +1)(\mu +2)]`$. Therefore, the small-$`v`$ asymptotics of the initial velocity distributions governs not only the scaling exponents but also the entire shape (including the large velocity tail) of the scaling function $`\mathrm{\Phi }(z)`$. This scaling behavior indicates that at time $`t`$, most cars moving initially with velocities larger than the typical velocity scale, $`v(t)t^\beta `$, have already joined clusters led by slower cars, while cars slower than this velocity scale are still driving with their initial velocity.

### B The size-velocity distribution

Given the nature of the model, a car is only affected by the initial configuration of cars ahead of it. This key feature enables solution of the cluster velocity distribution and it allows treatment of a more detailed quantity, the joint size-velocity distribution. To obtain $`P_m(v,t)`$, the density of clusters of size $`m`$ and velocity $`v`$, it is useful to introduce the cumulative distribution, $`Q_m(v,t)`$, the distribution of clusters of velocity $`v`$ containing at least $`m`$ cars. Knowledge of this cumulative distribution yields the joint size-velocity distribution via differencing, $`P_m(v,t)=Q_m(v,t)Q_{m+1}(v,t)`$.

Consider the first nontrivial quantity, $`Q_2(v,t)`$, the probability distribution of clusters of velocity $`v`$ with at least two cars at time $`t`$. This quantity is equal to the product of the probability that lead car has survived up to time $`t`$, $`P(v,t)`$, and the probability that the car trailing it actually experiences a collision prior to time $`t`$. Let $`x_1`$ and $`v_1`$ be the initial position and the initial velocity of the trailing car, respectively. For such a collision to occur, the trailing car must be faster than the lead car, $`v_1>v`$, and the interval separating the two cars must be initially free of other cars. The probability for this composite event is the product of the probabilities of each individual event. Given a random (Poisson) spatial distribution, the probability an interval is empty is exponential in its length, and the collision probability is

$$Q_2(v,t)=P(v,t)_v^{\mathrm{}}𝑑v_1P_0(v_1)_{x_1<(v_1v)t}𝑑x_1\mathrm{exp}(x_1).$$
(7)
The fact that the trailing car cannot be slowed down by any other car before colliding with the lead car is crucial in obtaining this solution. This solution can be generalized to arbitrary cluster sizes. Following the two-car cluster case, one simply integrates over all the initial positions and velocities of the consecutive cars to eventually collide with the lead car in the cluster to give

$$Q_m(v,t)=P(v,t)\underset{i=1}{\overset{m1}{}}_v^{\mathrm{}}𝑑v_iP_0(v_i)\underset{x_1+\mathrm{}+x_i<(v_iv)t}{}𝑑x_i\mathrm{exp}(x_i).$$
(8)
The integration limits reflect the fact that all the colliding cars must be faster than the lead car, and the restriction on the integration limits ensures that cars are sufficiently close to the lead car so that collisions indeed occur.

Given the cumulative car distribution, the joint size-velocity distribution can be formally obtained. Since $`mt^\alpha `$ and $`vt^\beta `$, we anticipate the following scaling behavior: $`P_m(v,t)t^{\beta 2\alpha }\mathrm{\Psi }(mt^\alpha ,vt^\beta )`$. This indeed holds and from Eq. (8) one obtains the scaled joint distribution

$$\mathrm{\Psi }(x,z)=cz^\mu (x+z)^{\mu +1}\mathrm{exp}\left[b(x+z)^{\mu +2}\right],$$
(9)
with $`c=a^2/(\mu +1)`$. Again only two parameters, $`a`$ and $`\mu `$ characterizing the small velocity characteristics of the initial conditions, are needed to fully describe the asymptotic state of the system. The joint distribution (9) provides a comprehensive description of the traffic clustering process. It may be considered as the counterpart of the well-known result for diffusion-controlled aggregation in one-dimension .

Integration of the scaling function with respect to $`x`$ reproduces the scaled velocity distribution $`\mathrm{\Phi }(z)`$. The complementary scaled size distribution cannot be found in a closed elementary form, except for the special case of asymptotically flat distributions ($`\mu =0`$) where both of the single variable scaling functions are purely Gaussian.

### C Generalizations

A natural generalization is to spatially heterogeneous initial velocity distributions, $`P_0(x,v)`$. The time and space dependent cluster velocity distribution, $`P(x,v,t)`$, follows from a straightforward generalization of the basic derivation in the homogeneous case

$$P(x,v,t)=P_0(xvt,v)\mathrm{exp}\left[_0^v𝑑v^{}_{xvt}^{xv^{}t}𝑑x^{}P_0(x^{},v^{})\right].$$
(10)
For instance, consider the special case where cars are uniformly distributed in the region $`x0`$ while the region ahead is empty. Here, one finds a governing length scale $`xvtt^\alpha `$ with the same exponent $`\alpha `$ as in Eq. (3). This length scale characterizes a propagating front of clusters, and the space dependent concentration $`c(x,t)`$ becomes a function of the scaling variable $`X=xt^\alpha `$, namely $`c(x,t)=t^\alpha 𝒞(X)`$. Far from the origin, the scaled density decays as $`𝒞(X)X^1`$ implying $`c(x,t)x^1`$ for $`xt^\alpha `$. Consequently, the total number of clusters in the originally empty region, $`N(t)=_0^{\mathrm{}}𝑑xc(x,t)`$, grows logarithmically slow with time

$$N(t)\mathrm{ln}t.$$
(11)
This growth law is universal as the dependence on the details of the initial velocity distribution is secondary, entering only via the prefactor.

In summary, a scaling asymptotic behavior characterizes the kinetics of clustering in no-passing zones of one lane roadways. The corresponding scaling exponents and scaling functions are characterized by the small-velocity statistics of the initial velocity distributions. Remarkably, it is possible to derive the exact Boltzmann equation in this case.

## III Passing zones

We now describe the complementary case of passing zones where fast cars can pass slow cars. The model we consider is a straightforward generalization of the no-passing case. The initial conditions are identical: cars are distributed randomly in space with concentration $`c_0`$ and their velocity is drawn from the intrinsic velocity distribution $`P_0(v)`$. The characteristic velocity scale is taken to be $`v_0`$. In the absence of other cars, cars drive ballistically with their intrinsic velocity. In the presence of other cars, two competing mechanism may cause a change in the car velocity. Collisions lead to slowing down: when a cluster overtakes a slower cluster, a larger cluster moving with the smaller of the two velocities forms. Passing leads to a velocity increase: every car inside a cluster may spontaneously pass the lead car and resume driving with its intrinsic velocity. The corresponding passing rate equals a constant, $`t_0^1`$. This is a significant simplification; in realistic situations only the first few trailing cars may be able to pass.

It proves convenient to introduce dimensionless velocity, space, and time variables: $`v/v_0v`$, $`xc_0x`$, $`c_0v_0tt`$. This rescales the passing rate, $`t_0^1R^1`$, where

$$R=\frac{t_{\mathrm{pas}}}{t_{\mathrm{col}}}=c_0v_0t_0$$
(12)
is the ratio of the passing time $`t_{\mathrm{pas}}=t_0`$ to the collision time $`t_{\mathrm{col}}=(c_0v_0)^1`$. We term this fundamental dimensionless quantity the “collision number” and denote it $`R`$ as it is reminiscent of the Reynolds number – the small $`R`$ limit is straightforward and the large $`R`$ limit is characterized by boundary layers.

The starting point for kinetic theory is again the cluster velocity distribution $`P(v,t)`$. The approximate Boltzmann equation reads

$`{\displaystyle \frac{P(v,t)}{t}}`$ $`=`$ $`R^1\left[P_0(v)P(v,t)\right]`$ (13)
$``$ $`P(v,t){\displaystyle _0^v}𝑑v^{}(vv^{})P(v^{},t).`$ (14)
This evolution equation assumes molecular chaos, namely that the stochastic passing events effectively mix the velocities, and therefore, spatial correlations can be neglected. This is clearly an approximation, as the collision integral does not coincide with the exact collision integral (2) derived in the no-passing case (the $`R\mathrm{}`$ limit). Still, this term reflects the fact that collisions occur only with slower clusters and that the collision rate is proportional to the velocity difference. The passing term is exact since the concentration of slowed down cars with intrinsic velocity $`v`$ equals $`P_0(v)P(v,t)`$.

In contrast with the no-passing case, the process is now stochastic in nature and the system approaches a nontrivial steady state. Setting the time derivative in Eq. (13) to zero we see that the steady state cluster velocity distribution $`P(v)P(v,t=\mathrm{})`$ satisfies the integral equation

$$P(v)\left[1+R_0^v𝑑v^{}(vv^{})P(v^{})\right]=P_0(v).$$
(15)
Given the intrinsic velocity distribution this relation gives the final cluster velocity distribution only implicitly. In contrast, the inverse problem is simpler as knowledge of the final distribution, the observed quantity in real traffic flows, gives explicitly the intrinsic distribution. We confirm that in the limit $`R\mathrm{}`$, all clusters move with the minimal velocity, while in the limit $`R0`$, all cars move with their intrinsic velocity $`P(v)P_0(v)`$.

The integral equation (15) can be transformed into a differential one using the auxiliary function $`Q(v)=R^1+_0^v𝑑v^{}(vv^{})P(v^{})`$, from which $`P(v)=Q^{\prime \prime }(v)`$. Thus Eq. (15) becomes

$$Q(v)Q^{\prime \prime }(v)=R^1P_0(v).$$
(16)
The boundary conditions are $`Q(0)=R^1`$ and $`Q^{}(0)=0`$. The auxiliary function $`Q(v)`$ gives a comprehensive description of the steady state. Calculation of important quantities such as the flux $`J`$ requires knowledge of $`G(v)`$, the car velocity distribution. This quantity satisfies $`1=𝑑vG(v)`$ and $`J=𝑑vvG(v)`$. Following an involved calculation that requires solution of a higher order velocity distribution , the car velocity distribution is derived explicitly in terms of $`Q(v)`$

$$G(v)=P(v)\left[1+R_v^{\mathrm{}}𝑑wP_0(w)_v^w\frac{du}{[RQ(u)]^2}\right].$$
(17)
Hence, for arbitrary intrinsic velocity distributions, the entire steady state problem is reduced to the nonlinear second order differential equation (16). Given $`Q(v)`$, steady state distributions such as $`P(v)`$ and $`G(v)`$ can be calculated using the explicit formulas above.

Except for a few special cases, one can not solve the differential equation (16) analytically. Nevertheless, the formal solution above can be used to evaluate generic features of the flow. The dimensionless collision number $`R`$ is extremely useful. For low collision numbers, a perturbation solution in powers of $`R`$ can be constructed, as the steady state differs weakly from the initial state. For high collision numbers, a boundary layer analysis is possible as sufficiently small velocities are not affected by collisions. These two limits are quantitatively analyzed as follows.

### A Low Collision Numbers

The flow characteristics in the collision-controlled regime, $`R1`$, can be analyzed systematically as a perturbation series in $`R`$. For example, the cluster velocity distribution and the car velocity distribution read

$`P(v)P_0(v)\left[1R{\displaystyle _0^v}𝑑v^{}(vv^{})P_0(v^{})\right],`$ (18)
$`G(v)P_0(v)\left[1+R{\displaystyle _0^{\mathrm{}}}𝑑v^{}(v^{}v)P_0(v^{})\right].`$ (19)
Consequently, average quantities such as the flux and the average cluster size vary linearly in $`R`$ in this free flow regime, $`J=J_0\mathrm{const}\times R`$, and $`m=1+\mathrm{const}\times R`$. The proportionality constant in the case of the flux equals the variance in the initial velocity distribution, indicating that the larger the initial velocity fluctuations, the larger the reduction in the flux.

Therefore, weakly interacting “laminar” flows arise in the $`R0`$ limit. Technically, the steady state remains close to the initial state and a perturbation series in the collision number is possible. Here, the assumptions made in our theory are justified, as the cluster sizes are small, and at the leading order, a simplified model were all cars in the cluster can pass coincides with a more realistic model where only the first few cars may pass. In fact, a basic prediction of the model, namely linear growth of the average cluster size with the flux is consistent with empirical data, obtained from observations of traffic flows in a secondary rural road in Los Alamos, New Mexico .

### B High Collision Numbers

The limit of high collision numbers corresponds to dense, congested flows where large clusters form. The analysis in this passing-controlled regime is more subtle since the condition $`R_0^v𝑑v^{}(vv^{})P_0(v^{})1`$ is satisfied only for small velocities. No matter how large $`R`$ is, sufficiently slow cars are not affected by collisions, and $`P(v)`$ is still given by Eq. (18) when $`vv^{}`$. The threshold velocity $`v^{}v^{}(R)`$ is estimated from $`R_0^v^{}𝑑v(v^{}v)P_0(v)1`$.

In the limit $`R\mathrm{}`$ limit, statistics of the slowest cars dominate the flow. Again, it is useful to consider intrinsic distributions with an algebraic small velocity form (5). For such distributions, the threshold velocity decreases with growing $`R`$ according to

$$v^{}R^{\frac{1}{\mu +2}}.$$
(20)
One can show that the flux is proportional to this velocity $`Jv^{}`$. For $`vv^{}`$, the collision integral in Eq. (15) dominates over the constant factor and $`RP(v)_0^v𝑑v^{}(vv^{})P(v^{})v^\mu `$. Anticipating an algebraic behavior for the cluster velocity distribution, $`P(v)R^\sigma v^\delta `$ when $`vv^{}`$, gives different answers dependent on whether the cutoff exponent $`\mu `$ is positive or negative. The leading behavior for $`vv^{}`$ can be summarized as follows

$$P(v)\{\begin{array}{cc}(v^{})^\mu (v/v^{})^{\mu 1}\hfill & \mu <0\text{;}\hfill \\ (v/v^{})^1[\mathrm{ln}(v/v^{})]^{\frac{1}{2}}\hfill & \mu =0\text{;}\hfill \\ (v^{})^\mu (v/v^{})^{\frac{\mu }{2}1}\hfill & \mu >0\text{.}\hfill \end{array}$$
(21)
On the other hand, $`P(v)P_0(v)`$ for $`vv^{}`$. This shows that the velocity distribution develops a boundary layer structure, the size of which vanishes in the infinite collision number limit. Inside the boundary layer, the velocity distribution is only marginally lower than its initial values, while the bulk of the velocities are strongly suppressed. Similar to the threshold velocity $`v^{}`$, macroscopic characteristics of the flow depend algebraically on $`R`$. For example, the average cluster size is

$$m\{\begin{array}{cc}R^{(\mu +1)/(\mu +2)}\hfill & \mu <0\text{;}\hfill \\ (R/\mathrm{ln}R)^{1/2}\hfill & \mu =0\text{;}\hfill \\ R^{1/2}\hfill & \mu >0\text{.}\hfill \end{array}$$
(22)
Two distinct regimes of behavior emerge. For $`\mu >0`$, car-cluster collisions dominate while for $`\mu <0`$ cluster-cluster collisions dominate. Interestingly, in the cluster-cluster dominated regime, $`mR^\alpha `$ with the scaling exponent $`\alpha =(\mu +1)/(\mu +2)`$ as in the no-passing case (3). Thus in the passing case the cutoff exponent $`\mu `$ also plays an important role in characterizing the behavior. Moreover, the steady state behavior is much richer than that found for the clustering kinetics.

Despite the simplifying assumptions, the model results in realistic behavior. The overall picture is both familiar and intuitive: due to the presence of slower cars, clusters form and the overall flux is reduced. For heavy traffic, the characteristics of the flow are solely determined by the distribution of slow cars. A single dimensionless parameter, the collision number $`R`$, ultimately determines the nature of the steady state.

## IV The Maxwell model

While a comprehensive analysis of the steady state velocity distributions is possible using the approximate kinetic theory (13), other important questions such as the relaxation toward the steady state and the nature of the cluster size distribution  remain unanswered. To address these issues we consider a further approximation where the collision rate is taken to be uniform . This approximation, known as the Maxwell model, is very useful in kinetic theory  and it has been recently applied to granular gases as well . In our case, it allows for a complete exact solution of the time dependent behavior, and additionally, it leads to closed evolution equations for the cluster-size distribution.

### A Relaxation

In the Maxwell approximation, the collision rate $`vv^{}`$ in the Boltzmann equation (13) is replaced by a constant factor which we set equal to unity. The corresponding rate equation for the cluster velocity distribution reads

$`{\displaystyle \frac{P(v,t)}{t}}`$ $`=`$ $`R^1\left[P_0(v)P(v,t)\right]`$ (23)
$``$ $`P(v,t){\displaystyle _0^v}𝑑v^{}P(v^{},t).`$ (24)
Again, the analysis is performed via a properly defined auxiliary function, $`Q(v,t)=_0^v𝑑v^{}P(v^{},t)`$. The constant collision rate results in simpler differential equations, that are only first order in the velocity. The analog of Eq. (16) is the integrable steady state equation $`Q(v)Q^{}(v)=R^1P_0(v)`$. The resulting steady state properties are governed by $`R`$, with a boundary layer structure in the large $`R`$ regime. The quantitative characteristics are somewhat different and for example the threshold velocity decays with $`R`$ according to $`v^{}R^{1/(\mu +1)}`$ rather than Eq. (20). If, however, the collision rate is properly chosen, namely set equal to $`v`$ rather than unity, we recover Eq. (20).

Furthermore, the complete time dependence can be obtained analytically by integrating the partial differential equation $`Q_t=R^1Q_vQQ_v`$. In general, the relaxation is exponential $`P(v,t)P(v,t=\mathrm{})f(v)e^{t/\tau (v)}`$, with $`\tau (v)=R\left[1+2RI_0(v)\right]^{1/2}`$ where $`I_0(v)=_0^v𝑑v^{}P_0(v^{})`$. The relaxation time depends on the velocity and the collision number according to

$$\tau (v)\{\begin{array}{cc}R\hfill & vv^{}\text{;}\hfill \\ [R/I_0(v)]^{1/2}\hfill & vv^{}\text{.}\hfill \end{array}$$
(25)
While small velocities are governed by practically fixed relaxation times, large velocities are characterized by velocity dependent decay rates. Furthermore, a large range of relaxation scales exists, $`R^{1/2}<\tau <R`$, with larger scales corresponding to smaller velocities. Further analysis shows that the same relaxation times underlie the car velocity distribution. We expect that while the predictions of the Maxwell model are only approximate, it correctly predicts the existence of a spectrum of relaxation time scales, and that the qualitative nature of the time dependent behavior generally holds.

### B The size distribution

The size distribution obeys closed evolution equations in the Maxwell model and can be solved exactly . It can also be used to address the nature of the passing mechanism. To demonstrate this we consider the model where only the next-to-leading car in the cluster may pass and resume driving with its intrinsic velocity. From numerical simulations of this model, we find two distinct phases. In the laminar regime, clusters are generally small, specifically the cluster size distribution is exponentially suppressed for sufficiently large sizes. When the collision number exceeds a certain threshold, an infinite cluster is formed, i.e., a finite fraction of the cars in the system are in the cluster behind the slowest car. Furthermore, in this jammed phase the size distribution of finite clusters has a fat tail close to a power-law, $`P_mm^\tau `$, with $`\tau 2`$.

In the Maxwell model framework, the cluster size distribution $`P_m(t)`$ obeys a closed system of rate equations

$`{\displaystyle \frac{dP_m}{dt}}=R^1[P_{m+1}P_m]cP_m+{\displaystyle \frac{1}{2}}{\displaystyle \underset{i+j=m}{}}P_iP_j,`$ (26)
$`{\displaystyle \frac{dP_1}{dt}}=R^1[P_2P_1+c]cP_1.`$ (27)
These equations were derived by enumerating all possible ways in which clusters evolve. For instance, consider Eq. (27). Collisions reduce the density of single cars, and the collision rate is clearly equal to $`c`$, as it is velocity-independent. The escape term in Eq. (27) is understood by observing that the rate of return of single cars into the system is $`2P_2+_{j3}P_j=P_2P_1+c`$. Here $`P_2`$ is singled out since passing transforms it into two single cars while an escape from larger clusters produces only one freely moving car.

Similar equations were previously studied in the context of aggregation-fragmentation processes . Utilizing the approach of Ref.  we find that a phase transition occurs at $`R_c=1`$ . For large $`m`$, the steady state size distribution is

$$P_m\{\begin{array}{cc}m^{3/2}\left[1(1R)^2\right]^m\hfill & R<1\text{,}\hfill \\ m^{5/2}\hfill & R1\text{.}\hfill \end{array}$$
(28)
Hence in the laminar regime, the size distribution decays exponentially in the large size limit. In the congested phase, the size distribution has a power law tail, and in addition there is an infinite cluster that contains the following finite fraction of cars in the system:

$$I=\{\begin{array}{cc}0,\hfill & R<1\text{;}\hfill \\ 1R^1,\hfill & R>1\text{.}\hfill \end{array}$$
(29)
Interestingly, this phase transition is similar to phase transitions in driven diffusive systems without passing . Furthermore, The formation of an infinite cluster is reminiscent of Bose-Einstein condensation .

## V Discussion

The most important question raised by the above results concerns the validity of the “mean-field” Boltzmann equation (13). Although passing is a stochastic mixing mechanism that diminishes correlations between the velocities and the positions of the cars, such correlations do exist, and it will be interesting to determine whether quantitative predictions such as the scaling behaviors (20) and (22) are altered by spatial correlations. Similarly, the collision term in Eq.(26) is written in a mean-field spirit and that may be the reason for the discrepancy between the theoretical prediction $`\tau =5/2`$ and the numerically observed value $`\tau 2`$ of the decay exponent $`P_mm^\tau `$.

The primary feature of our model is quenched disorder, which manifests itself in the random assignment of intrinsic velocities. Road conditions (construction zones, intermittent passing zones, turns, hills, etc.) present another source of quenched randomness in actual roads , which is ignored in our model. Quenched disorder significantly affects characteristics of many-particle systems, especially in low spatial dimensions . We have seen that this general conclusion clearly applies to our one-dimensional traffic model. Little is known analytically on the influence of the spatial disorder.

Finally, one may modify the passing rule so that when a car overtakes a slow car, it acquires a new velocity drawn from the distribution $`P_0(v)`$ rather than a pre-assigned velocity . This elementary zero-memory model remains highly non-trivial even in the collision-controlled limit $`R0`$ where clustering can be disregarded. The fate of the system is again determined by the behavior of the intrinsic velocity distribution near its lower cutoff. If $`P_0(v)`$ vanishes in this limit, the system reaches a steady state, otherwise, the system evolves indefinitely. Specifically, for intrinsic distributions with an algebraic small velocity tail (5) the long-time asymptotics of the average velocity reads

$$v(t)\{\begin{array}{cc}\mathrm{const}\hfill & \mu >0\text{;}\hfill \\ (\mathrm{ln}t)^1\hfill & \mu =0\text{;}\hfill \\ t^\mu \hfill & 1<\mu <0\text{.}\hfill \end{array}$$
(30)
These results were derived in a simplified Boltzmann framework. In particular, the most interesting behavior in the evolving regime was obtained by assuming that as $`t\mathrm{}`$, cars can be divided into two groups, the small group of “active” cars which move with velocities $`v1`$ and the vast majority of “creeping” cars that hardly move at all. We then ignored collisions between creeping cars (since their relative velocity is very small) and collisions between active cars (since their density is small). Thence, the velocity distribution of active cars obeys a linear Boltzmann-Lorentz equation which was solved to give (30). Comparison with results of molecular dynamics simulations suggests that the mean-field theory description is asymptotically exact. It will be interesting to confirm this result rigorously.

We are thankful to our collaborators Slava Ispolatov and Sid Redner. This research was supported by DOE (W-7405-ENG-36) and NSF(DMR9978902).


text:
# The Spectrum of Earthshine: A Pale Blue Dot Observed from the Ground

## 1 Introduction

Earthshine is sunlight which has reflected from the Earth onto the dark side of the moon and back again to Earth. Here “dark side” refers to that portion of the lunar surface which, at any instant, faces the Earth but does not face the sun. Both earthshine, from the dark side of the moon, and moonlight, from the bright side of the moon, are transmitted through the same air mass just prior to detection, and thus suffer the same extinction and imposed absorption features. Their ratio is the reflection coefficient of the Earth multiplied by a geometric factor, and by the ratio of phase effects of the moon. (We ignore weak additive Raman-shifted spectra from both bodies, i.e., the Ring effect.) We report here the results of our effort to obtain a spectrum of the Earth as it would be seen from outside the solar system. Our goal is to pave the way for interpreting extrasolar terrestrial planet spectra which might be obtained in the future, for example as part of the Terrestrial Planet Finder (TPF) program.

## 2 Observations

We observed the spectrum of earthshine and moonlight soon after sunset on June 24/25 2001 (25/26 UT) from Kitt Peak, Arizona, using the Steward Observatory 2.3m telescope with its Boller and Chivens spectrograph, a CCD detector and a 300 line/mm grating. The raw spectrum is shown in Fig. 1. A spectral range from 4800 to 9200Å was recorded, at a resolution of about 600, or a spectral width of about 8-9Å. The telescope drive was set to lunar rate and was guided in right ascension, keeping the limb in the same position on the slit to within about a minute of arc in declination between the different observations. At the time of these observations, about 1/4 of the moon appeared sunlit, and 3/4 was nominally dark.

The 1.5 arcsec wide spectrograph slit has a length of about 200 arcsec. It was set perpendicular to the lunar limb, with approximately half of its length on the limb and half on adjacent sky. Both the sunlit and Earth-lit limb spectra were corrected for scattered light in the telescope by subtracting the adjacent sky spectrum; this amounted to about 2/3 of the total signal in the earthshine spectrum. Each earthshine spectrum was bracketed by moonlight spectra.

Both nights were partly cloudy, and the observations were made during breaks in the clouds. Data from the first night have more noise than data from the second night, but are in qualitative agreement. The second night produced two independent sets of observations; the first was taken under more clear conditions, and is the one reported here.

The measured reflectivity of the integrated Earth is shown in Fig. 2. The reflectivity is the ratio of the earthshine to moonlight spectra, both sky-subtracted. We have also corrected for the color effect of lunar phase (Lane and Irvine, 1973), which makes the apparent reflectivity about 16% too strong in the blue; this color bias was removed by multiplying the apparent reflectivity by a 8-segment piecewise-linear factor which ranged from 0.835 at 5000Å to 1.00 at 9000Å.

The reflectivity in Fig. 2 has also been corrected for CCD fringing. Dome flats were not adequate to exactly match the lunar-spectrum fringing in the thinned CCD chip because flexure in the instrument and/or telescope changed the position of the fringes slightly. Use of a quartz lamp internal to the spectrograph ensured that fringes were at the same location on the chip as for the lunar observations, but their amplitude was suppressed because the lamp does not illuminate the slit with an f/9 beam as does light that passes through the telescope. We used the internal lamp to partially correct the fringes, by dividing by the lamp spectrum, but a noticeable fringe pattern remained. Five residual interference fringe packets in the red part of the spectrum were least-squares fitted with smooth gaussian envelope sinusoids and subtracted from the data. The fitting procedure removed most of the remaining fringing, but there remains a residual noise which is larger in this region than elsewhere in the spectrum. The zero-mean subtracted component is explicitly displayed (offset in intensity) beneath the data in Fig. 2.

## 3 Spectrum Model and Analysis

Our model of the Earth’s reflection spectrum is a simple box model in which each box represents a potentially major contributor to the final spectrum. By using a box model, we ignore the interactions between components (e.g., reflections between clouds and ground, multiple scattering, etc.); in other words, we focus on the net resulting spectrum, not the interplay of generating mechanisms. The model has 9 independent parameters, each of which scales an independent physical component. The standard-atmosphere component is discussed in Traub and Jucks (2002), and the other components are discussed briefly here.

(1) Optical depth of standard atmospheric species H2O, O2, and O3. The optical depth is calculated at very high spectral resolution, using a model atmosphere with thin layers from 0 to 100 km altitude, standard mixing ratio profiles of species, and Doppler and pressure broadening. The resulting spectrum is smeared and sampled to match the experimental wavelength grid and resolution.

(2) Optical depth of stratospheric O3. Here O3 is the only absorber, and atmospheric layers are limited to stratospheric altitudes, otherwise similar to the full standard atmosphere component.

(3) Rayleigh scattering. This is a source term with $`\lambda ^4`$ functional form.

(4) Aerosol scattering. This is also a source term but with $`\lambda ^{1.3}`$ form.

(5) Reflectivity of high clouds. This is a constant reflectivity term, independent of wavelength, as from a very high cloud.

(6) Reflectivity of ocean water. This is light which has penetrated the surface and has been scattered out of the water. It has a $`10^{(\lambda \lambda _0)/\mathrm{\Delta }\lambda }`$ form, giving the ocean its blue color, where $`\lambda _0=5500`$Å and $`\mathrm{\Delta }\lambda =1300`$Å here.

(7) Reflectivity of ocean pigment. This is ocean pigment due to chlorophyll and its products, here represented by a gaussian centered at 5500Å and with full width at half maximum 650Å, giving waters with phytoplankton a green color. The longer-wavelength chlorophyll feature (item 8) is not seen in practice; water absorbs quite strongly at these wavelengths (item 6).

(8) Reflectivity of vegetated land. This is the reflectivity of vegetation containing chlorophyll, on land, dominated by a sharp rise in reflectivity for wavelengths longer than about 7200Å, plus some smaller bumps at shorter and longer wavelengths. It is modeled here by functional forms approximating these features. (If our eyes were sensitive to wavelengths longer than 7200Å, all vegetation would appear to be infrared bright, not green!)

(9) Reflectivity of low clouds and surface, below bulk of atmosphere. The nominally featureless reflectivity of land, low clouds, and the 2% Fresnel reflectivity of the air-ocean interface are lumped into this single term, independent of wavelength.

In terms of our model, the two dominant components in the spectrum are the “high cloud” continuum and the “clear atmosphere” spectrum of water, ozone, and oxygen. All components are displayed in Fig. 2, where the linear sum of the 7 component spectra is equal to the smooth model function superposed on the data. The main features of water, ozone, and oxygen are labeled; weak bands of water are not labeled. The relative strength of each spectral component is proportional to the intrinsic reflectivity times the effective projected area on Earth of that component. For example, suppose that the clear atmosphere spectrum arises entirely over water, which has a Fresnel reflectivity of 0.02, and that the continuum spectrum is from clouds, with a nominal reflectivity of 0.30. Then the fact that we observe these components to have similar strengths indicates that the clear atmosphere area is roughly 15 times greater than the cloudy area, i.e., the model Earth has a relatively clear sky.

Rayleigh scattered light from molecules is the next largest component. In our box model this component is multiplied by the transmission of the clear atmosphere, which is how it is displayed in Fig. 2. The remaining components are relatively small, and carry less significance than the dominant three.

The vegetation “red edge” signal in Fig. 2 is the well-known signature of chlorophyll in land vegetation, described above. The earthshine spectrum shows increased reflectivity in the 7200 to 7900Å region, suggestive of a vegetation signal in this range. However the data do not show continued high reflectivity at longer wavelengths (7900 to 9200Å), as would be expected for a vegetation reflection profile, making the interpretation ambiguous. Nevertheless, on the basis of the abrupt reflectivity increase at 7200Å we suspect that the vegetation signal is indeed present, and that either the data or the model fail at longer wavelengths, masking this signal.

The three remaining components are weakest. The ocean clear-water signal, the aerosol signal, and the ocean phytoplankton signal all make relatively small contributions to the spectrum. Note that there is a superficial similarity between each of the blue continuum spectra (Rayleigh, ocean, and aerosol), however their wavelength variations are quite different when viewed over a wide range of wavelengths, and this allows the components to be independently determined, with relatively small correlation.

In summary, the earthshine spectrum shows strong evidence for cloud, ozone, water, oxygen, Rayleigh scattering, and very probably vegetation signals. Apparently weak components are aerosols, ocean water, and phytoplankton signals.

## 4 Discussion

In order to characterize an extrasolar terrestrial planet using visible light, there are three main avenues of approach. (1) The diurnal photometric variability could be measured. Ford et al. (2001) have shown that for an Earth-like planet, variation on the order of 100% may be expected. (2) The spectral reflectivity of the atmosphere could be measured. As shown here and in Traub and Jucks (2002), absorption features due to ozone, water vapor, and molecular oxygen could be seen, as well as enhanced reflectivity due to Rayleigh and aerosol scattering. (3) The spectral reflectivity of the surface could be measured. As we have shown in Fig. 2. This potentially includes the reflectivity of clear ocean water (blue), phytoplankton in the ocean (green), and land vegetation (red).

The vegetation signal in Fig. 2 is relatively small, about 6% above the nearby continuum. But at the time of observation only about 17% of the projected area was land, as shown in the simulated cloud-free Earth image111Earth simulation images are available at http://www.fourmilab.ch/earthview/vplanet.html in Fig. 3. Of this fraction of land, only part is cloud-free and vegetated, so in principle the vegetation signal could be at least 6 times larger than we observe. An infrared view of the whole earth222http://www.ssec.wisc.edu/data/comp/ir/ showed little or no cloud cover over the land masses illustrated in Fig. 3, so it is reasonable to expect that we might be seeing a large part of the available vegetation signal in our spectrum.

For example, this expectation is borne out in data taken by the GOME satellite333For a sample spectrum, see http://earth.esa.int/symposia/data/popp in which the spectral region from about 3100 to 8000Å was observed in a nadir-pointing mode, and the field of view was relatively small, and filled with vegetated land in northern Italy. In that reflection spectrum, the vegetation signal is about 2.4 times larger than the adjacent continuum, showing that this is a potentially very large differential signal on a similarly vegetated extrasolar planet.

The residual spectral mismatch in Fig. 2 in the range 7300-7900Å suggests that yet more vegetation signal could be accommodated in the model, but that this may be thwarted by a bias in either the data or the model in the 7900 to 9200Å region. Although a more detailed reflectivity model of the composite Earth has yet to be developed, we believe that the enhancement of reflectivity redward of 7300Å is significant, and is an indicator of vegetation.

## 5 Comments

Our spectrum of the Earth as it would appear to an extrasolar observer will be useful for learning how to analyze the spectra of extrasolar planets, discovering the key spectral features of Earth analogs, and determining the necessary spectral resolution and signal-to-noise ratios. This spectrum illustrates both atmospheric and surface reflectivity features. It was not feasible to look for variability because of the limited observing period.

We note that the reflectivity of the Earth has not been static throughout the past 4.5 Gyr. In particular, oxygen and ozone become abundant perhaps only 2.3 Gyr ago, affecting the atmospheric absorption component of the reflection spectrum. Then perhaps about 2.0 Gyr ago a green phytoplankton signal developed in the oceans. About 0.44 Gyr ago, an extensive land plant cover developed, generating the red chlorophyll edge in the reflection spectrum (Lunine, 1999).

An observer in a nearby stellar system would have been able to use the oxygen and ozone absorption features to suspect the presence of life on Earth anytime during the past 50% of the age of the solar system. The extrasolar observer would have been able to use the chlorophyll red-edge reflection feature to confirm the presence of life on Earth anytime during the most recent 10% of the age of the solar system.

As the Earth evolved there have also been significant changes in the temperature structure, the abundances of H2O, CO2, CH4, and perhaps the character of cloud coverage. Some of these changes may have been sudden and even oscillatory, as opposed to gradual. Thus the spectrum of the evolving Earth has shown dramatic changes (e.g., Traub and Jucks (2002)), and we may expect that similar signatures will be found on extrasolar terrestrial planets.

We suggest that further observations of Earthshine should be made. One purpose is to validate our ability to characterize extrasolar terrestrial planets using visible wavelengths In particular, it would be valuable to make earthshine observations from space as well as ground, to remove the complicating factor of observing through the Earth’s atmosphere. It would also be valuable to observe different types of land and ocean and cloud cover on Earth, and furthermore to do so for a full diurnal cycle in order to test the spectroscopy and the prediction of variable photometry.

We thank our colleagues on the TPF biomarkers team study (Des Marais et al., 2002) initiated by Charles Beichman and led by David Des Marais. We thank Felix Kogan, James Marr, Sara Seager, Jim Kasting, and an anonymous referee for advice. This work was supported in part by funds from Jet Propulsion Laboratory’s TPF studies, to NW via the Lockheed-Martin corporation, and to WT and KJ via the Ball Aerospace corporation.


text:
# Stationary cantilever vibrations in the oscillating cantilever-driven adiabatic reversals – magnetic resonance force microscopy technique

## I Introduction

Magnetic resonance force microscopy (MRFM) based on a cyclic adiabatic inversion (CAI) is considered as one of the most promising roads to the ultimate goal of a single-spin detection in solids (see, for example, 2 ; 3 ). Typically CAI is generated using the frequency modulation of the external radio-frequency (rf) field. In this case, a paramagnetic moment of a sample follows the effective magnetic field in the rotating system of coordinates (RSC), and influence the cantilever vibrations.

Recently a new technique called “oscillating cantilever-driven adiabatic reversals” (OSCAR) has been suggested and implemented in 1 . In this technique, the cantilever driven by an external force causes the CAI of the paramagnetic moment of a sample. The back reaction of the paramagnetic moment causes the frequency shift of the cantilever vibrations, which is supposed to be detected. The main purpose of this paper is a theoretical analysis of the stationary vibrations of the cantilever in the OSCAR technique. Our consideration is based on the classical equations of motion for the spin-cantilever system.

The paper is organized as follows. In Sec. II we introduce the model. The linear OSCAR regime is considered in Sec. III, and nonlinear regime is analyzed in Sec. IV. A perturbative approach and numerical results are presented in Sec. V. In Sec. VI, we analyze the damped oscillations of the cantilever in the absence of the external force. In Sec. VII we give a brief summary of our results.

## II Hamiltonian and equations of motion

A schematic setup of the studied system is shown in Fig. 1. A spherical ferromagnetic particle with magnetic moment, $`\stackrel{}{m}_F`$, is attached to the cantilever tip. A small paramagnetic cluster with magnetic moment, $`\stackrel{}{\mu }`$, which must be detected, is placed on the surface of non-magnetic sample beneath the tip of the cantilever. The whole system is placed into the high permanent magnetic field, $`\stackrel{}{B}_0+\mathrm{\Delta }\stackrel{}{B}`$, oriented in the positive $`z`$-direction. The external force, $`F(t)`$, drives the cantilever vibrations along the $`z`$-axis. The transversal rotating magnetic field, $`\stackrel{}{B}_1(t)`$, is applied to the paramagnetic cluster. We place the origin of our coordinate system at the equilibrium position of the cantilever tip.

We consider the cantilever tip as an oscillator with the effective mass, $`m^{}`$, and the effective spring constant, $`k_s`$. The classical Hamiltonian for the cantilever with the ferromagnetic particle and the paramagnetic cluster has the form,

$$=\frac{p_z^2}{2m^{}}+\frac{k_sz^2}{2}zF(t)\frac{\mu _0m_F}{2\pi (d+z)^3}\mu _z\stackrel{}{\mu }(\stackrel{}{B}_0+\mathrm{\Delta }\stackrel{}{B}+\stackrel{}{B}_1),$$
(1)
where $`p_z`$ and $`z`$ are the momentum and coordinate of the cantilever tip, $`\mu _0`$ is the permeability of the free space. Putting $`F(t)=F_0\mathrm{cos}(\nu t)`$ and taking into consideration the finite quality factor, $`Q`$, of the cantilever we write the equation of motion for the cantilever,

$$\ddot{z}+\omega _c^2z+\frac{q\mu _z}{(d+z)^4}+\frac{\omega _c}{Q}\dot{z}=f_0\mathrm{cos}(\nu t),$$
(2)
where $`\omega _c=(k_s/m^{})^{1/2}`$ is the unperturbed cantilever frequency, $`f_0=F_0/m^{}`$, and,

$$q=\frac{3\mu _0m_F}{2\pi m^{}}.$$
(3)
Next, we assume that the rf field $`\stackrel{}{B}_1`$ rotates in the $`(x,y)`$ plane with the frequency,

$$\omega _0=\gamma [B_0+B_d(0)].$$
(4)
Here $`\gamma `$ is the gyromagnetic ratio of the paramagnetic cluster, $`B_d(z)`$ is the dipole magnetic field produced by the ferromagnetic particle at the point of location of the paramagnetic cluster,

$$B_d(z)=\frac{\mu _0m_F}{2\pi (d+z)^3},$$
(5)
and $`B_d(0)`$ is the value of $`B_d(z)`$ at the equilibrium position of the cantilever, $`z=0`$. The equation of motion for the paramagnetic moment $`\stackrel{}{\mu }`$ in the RSC has the form,

$$\dot{\stackrel{}{\mu }}=\gamma [\stackrel{}{\mu }\times \stackrel{}{B}_{eff}].$$
(6)
Here $`\stackrel{}{B}_{eff}`$ is the effective magnetic field in the RSC with the $`x`$-component $`B_1`$ and the $`z`$-component $`\mathrm{\Delta }B+B_d^{}(z)`$, where $`B_d^{}(z)`$ is the oscillatory part of the dipole field produced by the ferromagnetic particle on the cluster:

$$B_d^{}(z)=B_d(z)B_d(0).$$
(7)
## III The linear OSCAR regime: small oscillations of $`\stackrel{}{\mu }`$

In this section we consider the linear OSCAR regime. Suppose that initially an auxiliary $`\pi /2`$-pulse changes the direction of the paramagnetic moment, $`\stackrel{}{\mu }`$, from $`+z`$ to $`+x`$ of the RSC. We also assume that the oscillatory part of the dipole field, $`B_d^{}(z)`$, is small compared to the rf field, $`B_1`$. Certainly we assume that the unperturbed cantilever frequency, $`\omega _c\gamma B_1`$, to keep the conditions of CAI. In quasi-static approximation a paramagnetic moment, $`\stackrel{}{\mu }`$, follows the effective field $`B_{eff}`$. Putting in (6) $`\dot{\stackrel{}{\mu }}=0`$, we obtain for $`|z|d`$:

$$\mu _x(t)\mu ,\mu _y(t)=0,\mu _z(t)=\frac{\mu }{B_1}\left[\mathrm{\Delta }B\frac{3\mu _0m_Fz(t)}{2\pi d^4}\right].$$
(8)
These equations describe small (linear) oscillations of $`\stackrel{}{\mu }`$ near the $`x`$-axis. Substituting the last expression in Eq. (8) into Eq. (2) we derive an approximate equation for the cantilever oscillations,

$$\ddot{z}+\omega _c^2z+\frac{\omega _c^{}}{Q^{}}\dot{z}=f_0\mathrm{cos}(\nu t).$$
(9)
Here,

$$\omega _c^{}=\omega _c+\mathrm{\Delta }\omega _c,\mathrm{\Delta }\omega _c=\frac{3\mu _0m_F\mu }{\pi m^{}\omega _cB_1d^5}\left(\mathrm{\Delta }B+\frac{3\mu _0m_F}{8\pi d^3}\right),Q^{}=Q\left(1+\frac{\mathrm{\Delta }\omega _c}{\omega _c}\right).$$
(10)
Equation (9) describes the motion of the linear oscillator with the effective frequency, $`\omega _c^{}`$, and the effective quality factor, $`Q^{}`$. Due to the back reaction of the paramagnetic moment on the cantilever the effective frequency and the quality factor of the cantilever depend on the permanent magnetic field, $`\mathrm{\Delta }B`$, (in our approximation, $`\mathrm{\Delta }BB_1`$). If $`\mathrm{\Delta }B>3\mu _0m_F/(8\pi d^3)`$, then both the frequency and the quality factor of the cantilever decrease. In the opposite case they increase.

## IV Nonlinear adiabatic regime: adiabatic reversals of $`\stackrel{}{\mu }`$

To increase the back reaction of $`\stackrel{}{\mu }`$ it is important to provide large oscillations (adiabatic reversals) of the paramagnetic moment. In this section we consider stationary vibrations of the cantilever in the nonlinear OSCAR regime. It is convenient to write the equations of motion in the dimensionless form:

$$Z^{\prime \prime }+Z+\frac{\lambda M_z}{(1+\alpha Z)^4}+\frac{1}{Q}Z^{}=\frac{1}{Q}\mathrm{cos}[(1+\rho )\tau +\vartheta _0],$$

$$M_{x}^{}{}_{}{}^{}=(\delta \chi Z)M_y,$$
(11)
$$M_{y}^{}{}_{}{}^{}=\epsilon M_z(\delta \chi Z)M_x,$$

$$M_{z}^{}{}_{}{}^{}=\epsilon M_y,$$

where we introduced the dimensionless time $`\tau =\omega _ct`$; prime means a differentiation over $`\tau `$, $`Z=z/A`$ is the dimensionless coordinate, $`A=f_0Q/\omega _c^2`$ is the unperturbed (in the absence of the magnetic moment $`\stackrel{}{M}`$) amplitude of the stationary cantilever vibrations in the resonant regime (when $`\omega =\omega _c`$), $`\stackrel{}{M}=\stackrel{}{\mu }/\mu `$ is the dimensionless magnetic moment, $`\delta =\gamma \mathrm{\Delta }B/\omega _c`$. The parameter $`\alpha =A/d`$ is small, $`\alpha 0.01`$. The dynamics is controlled by the following dimensionless parameters:

$$\lambda =\frac{3\mu _0m_F\mu }{2\pi d^4QF_0},$$

$$\chi =\frac{3\gamma \mu _0m_FQf_0}{2\pi \omega _c^3d^4},$$
(12)
$$\epsilon =\frac{\gamma B_1}{\omega _c},\rho =\nu /\omega _c1.$$

Suppose that the paramagnetic moment, $`\stackrel{}{M}`$, points initially in the direction of the effective magnetic field, $`\stackrel{}{B}_{eff}`$, and the cantilever points in the opposite direction, $`Z(0)=1`$. In this case the quasi-static motion of $`\stackrel{}{M}`$ is given by the expressions,

$$M_x(\tau )=\frac{\epsilon }{\sqrt{\epsilon ^2+(\delta \chi Z)^2}},$$

$$M_y=0,$$
(13)
$$M_z(\tau )=\frac{\delta \chi Z}{\sqrt{\epsilon ^2+(\delta \chi Z)^2}}.$$

Substituting (13) into the first equation in (11) we obtain the nonlinear equation for $`Z`$:

$$Z^{\prime \prime }+Z\frac{\lambda \chi Z}{\sqrt{\epsilon ^2+(\chi Z)^2}}+\frac{1}{Q}Z^{}=\frac{1}{Q}\mathrm{cos}[(1+\rho )\tau +\vartheta _0],$$
(14)
where we neglected the term $`\alpha Z`$ in the denominator in the third term in the left-hand side and put $`\delta =0`$. The third term in Eq. (14) corresponds to the modification of the potential energy of the cantilever, due to the interaction with the magnetic moment, $`\stackrel{}{M}`$, by the value,

$$\delta U(Z)=\frac{\lambda }{2\chi }\sqrt{\epsilon ^2+(\chi Z)^2}.$$

Now, we present an approximate “semi-quantitative” analysis of the stationary oscillations described by Eq. (14). The solution for the stationary driven oscillations of the cantilever, described by Eq. (14), can be written in the form,

$$Z=a(\rho )\mathrm{sin}[(1+\rho )\tau +\vartheta _0].$$
(15)
We define the frequency shift which corresponds to the shift of the maximum, $`a_{max}(\rho _1)`$, of the amplitude, $`a=a(\rho )`$, caused by the ferromagnetic sample. In order to estimate $`a_{max}`$ we replace,

$$Z^2\mathrm{sin}^2[(1+\rho )\tau +\vartheta _0]=\frac{1}{2}\{1\mathrm{cos}[2(1+\rho )\tau +2\vartheta _0]\}\frac{1}{2},$$

in the denominator in the third term in Eq. (14), and neglect the term $`\mathrm{cos}[2(1+\rho )\tau +2\vartheta _0]`$ because it is non-resonant. Then Eq. (14) takes the form,

$$Z^{\prime \prime }+\left[1\frac{\lambda \chi }{\sqrt{\epsilon ^2+\frac{\chi ^2}{2}}}\right]Z+\frac{1}{Q}Z^{}=\frac{1}{Q}\mathrm{cos}[(1+\rho )\tau +\vartheta _0].$$
(16)
The position, $`\rho _1`$, of the maximum of the amplitude, $`a_{max}`$, of the driven oscillations (the frequency shift) is,

$$\rho _1\rho _0=\frac{\lambda \chi }{2\sqrt{\epsilon ^2+\frac{\chi ^2}{2}}}\frac{\lambda }{\sqrt{2}},$$
(17)
where $`\rho _0=1/4Q^2`$ is the position of the maximum of the amplitude in the absence of the paramagnetic sample, and we suppose that $`\chi \epsilon `$.

For estimation of the value of the frequency shift the following parameters where used: $`D=1.5\times 10^7`$m is the diameter of the ferromagnetic particle with the volume $`V=1.8\times 10^{21}`$ m3, $`\mu _0m_F/V1.1`$T, $`k_c10^3`$N/m, $`\omega _c/2\pi 10^5`$Hz, $`A1`$nm, $`d100`$nm, $`B_110^3`$T. For these values of parameters we obtain,

$$\epsilon 280,\chi 2.5\times 10^3,\lambda 8.5\times 10^5(\mu /\mu _B),\alpha =0.01,$$
(18)
where $`\mu /\mu _B`$ is the paramagnetic moment expressed in units of the Bohr magneton. The corresponding frequency shift is,

$$\rho _1\rho _06\times 10^5(\mu /\mu _B).$$
(19)
This gives the frequency shift -6 Hz per one Bohr magneton.

## V Perturbation approach

The qualitative estimation presented above can be supported by application of the approach based on the perturbation theory developed by Bogoliubov and Mitropolskii in 4 . We look for the solution of Eq. (14) in the form,

$$Z=a(\tau )\mathrm{cos}[\psi ]+\lambda u_1(a,\psi ),$$
(20)
where $`\psi =(1+\rho )\tau +\vartheta (\tau )`$. The function $`u_1(a,\psi )`$ is the sum of the Fourier terms with the phases $`3\psi `$, $`5\psi `$, $`7\psi `$, $`\mathrm{}`$. The amplitudes of these terms decrease with increasing the Fourier number, $`n`$, as $`1/(2n+1)^2`$. The first non-vanishing term is small and equals to $`u_1(a,\psi )0.02\mathrm{cos}(3\psi )`$. This allows us to neglect the contribution of $`u_1(a,\psi )`$ into the expression for $`Z`$ in Eq. (20).

The slow varying amplitude, $`a(\tau )`$, and the phase, $`\vartheta (\tau )`$, in the first order of the perturbation theory satisfy the two coupled differential equations,

$$\frac{da}{d\tau }=\frac{\lambda }{2\pi }_0^{2\pi }\frac{\chi a\mathrm{cos}\psi \mathrm{sin}\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}\frac{a}{2Q}\frac{1}{Q(2+\rho )}\mathrm{sin}\vartheta ,$$
(21)
$$\frac{d\vartheta }{d\tau }=\frac{1}{8Q^2}\rho \frac{\lambda }{2\pi a}_0^{2\pi }\frac{\chi a\mathrm{cos}^2\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}\frac{1}{aQ(2+\rho )}\mathrm{cos}\vartheta .$$
(22)
Note, that the integral in the right-hand side of Eq. (21) is equal to zero. The integral in the right-hand side of Eq. (22) can be expressed through the elliptic integrals as (see Appendix),

$$4_0^{\pi /2}\frac{\chi a\mathrm{cos}^2\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}=4\left[\frac{1}{k}E(k)p^2kK(k)\right],$$
(23)
where $`k=1/\sqrt{1+p^2}`$, $`K(k)`$ and $`E(k)`$ are the complete elliptic integrals, respectively, of the first and second kind, $`p=\epsilon /(a\chi )`$. When $`p^21`$ one can decompose $`K(k)`$ and $`E(k)`$ as,

$$K(k)C+(C1)\frac{k_{}^{}{}_{}{}^{2}}{4}+\mathrm{},E(k)1+\left(C\frac{1}{2}\right)\frac{k_{}^{}{}_{}{}^{2}}{2}+\mathrm{},$$
(24)
where $`k_{}^{}{}_{}{}^{2}=1k^2p`$, $`C=\mathrm{ln}(4/k^{})\mathrm{ln}(4/p)`$. From Eqs. (23) and (24) we find the value of the integral in Eq. (22) for $`p1`$,

$$\frac{\lambda }{2\pi a}_0^{2\pi }\frac{\chi a\mathrm{cos}^2\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}\frac{2\lambda }{\pi a}\left[1\frac{p^2}{4}\left(2\mathrm{ln}\frac{4}{p}1\right)\right].$$
(25)
Substituting Eq. (25) to Eq. (22) we obtain,

$$\frac{da}{d\tau }=\frac{a}{2Q}\frac{1}{Q(2+\rho )}\mathrm{sin}\vartheta ,$$
(26)
$$\frac{d\vartheta }{d\tau }=\frac{1}{8Q^2}\rho \frac{2\lambda }{\pi a}\left[1\frac{p^2}{4}\left(2\mathrm{ln}\frac{4}{p}1\right)\right]\frac{1}{aQ(2+\rho )}\mathrm{cos}\vartheta .$$

We now calculate the position of the maximum of the amplitude, $`a(\rho )`$, in the stationary regime of driven oscillations using Eq. (26), and compare it with the results obtained in Sec. III. In the regime of driven oscillations $`a=const`$, $`\vartheta =const`$, and we must solve the system of two equations (26), where $`da/d\tau =0`$ and $`d\vartheta /d\tau =0`$. Canceling the phase, $`\vartheta `$, we have,

$$\frac{1}{a^2(2+\rho )^2}=\frac{1}{4}+Q^2\left(\frac{1}{8Q^2}+\rho +\frac{2\lambda }{\pi a}\right)^2,$$
(27)
where we neglected the term proportional to $`p^21`$. The amplitude, $`a`$, can be written as $`a=1+\beta `$, where $`\beta 1`$, so that,

$$\frac{1}{a(2+\rho )}=\frac{1}{(1+\beta )(2+\rho )}\frac{1}{2}\left(1\beta \frac{\rho }{2}\right).$$
(28)
Taking the square root from the both sides of Eq. (27) and using Eq. (28) we obtain,

$$\beta \frac{\rho }{2}2Q^2\left(\frac{1}{8Q^2}+\rho +\frac{2\lambda }{\pi }\right)^2,$$
(29)
where we put $`a1`$ in the denominator of the term proportional to $`\lambda `$ (i.e. we neglected the term of the order of $`\beta \lambda `$). The maximum of the function, $`\beta =\beta (\rho )`$, can be found from the condition $`d\beta (\rho _1)/d\rho =0`$ which yields,

$$\rho _1=\frac{1}{4Q^2}\frac{2\lambda }{\pi }.$$
(30)
This is approximately the same value as that given by Eq. (18), obtained from the qualitative considerations. The second term in Eq. (30) describes the influence of the paramagnetic moment reversals on the resonance frequency of the cantilever.

To verify our analytical results we solved numerically the exact equations of motion (11). Fig. 2 (solid line) demonstrates the dependence of the stationary amplitude of the cantilever vibrations, $`a`$, on the frequency detuning, $`\rho `$. (The stationary amplitude is achieved at $`\tau Q`$.) The initial conditions are taken in the form,

$$Z(0)=1,\dot{Z}(0)=0,M_x(0)=\frac{\epsilon }{\sqrt{\epsilon ^2+\chi ^2}},M_y(0)=0,M_z(0)=\frac{\chi }{\sqrt{\epsilon ^2+\chi ^2}}.$$
(31)
For these initial conditions at $`\tau =0`$ the paramagnetic moment, $`\stackrel{}{M}`$, points in the direction of the effective magnetic field, while the cantilever is displaced in $`z`$-direction ($`\vartheta _0=3\pi /2`$) from its equilibrium position.

Fig. 3 demonstrates the motion of the paramagnetic moment, $`\stackrel{}{M}(\tau )`$. One can see the close correspondence between the analytical and numerical solutions. Note that the frequency shift caused by the adiabatic reversals changes its sign if the paramagnetic moment points initially in the direction opposite to the effective magnetic field (while $`Z(0)=1`$). Dotted line in Fig. 2 depicts this case. We also should note, that decreasing the parameter $`\epsilon `$ (the $`x`$-component of the effective magnetic field) leads to the violation of the CAI conditions. Fig. 4 demonstrates this situation for $`\epsilon =28`$.

## VI Damped oscillations of the cantilever

The influence of the sample on the cantilever can be measured if one turns off the external force acting on the cantilever, and measures the frequency of small damped oscillations of the cantilever. In absence of the paramagnetic moment $`\stackrel{}{M}`$ the frequency of the oscillations is independent of time and equals to $`\sqrt{11/(4Q^2)}`$.

We look for the solution of the cantilever vibrations in the form: $`Z=a(\tau )\mathrm{cos}[\tau +\vartheta (\tau )]`$. Then the dynamical equations for the slow varying amplitude, $`a(\tau )`$, and phase, $`\vartheta (\tau )`$, in the presence of the sample and in the absence of the external force take the form,

$$\frac{da}{d\tau }=\frac{a}{2Q},$$
(32)
$$\frac{d\vartheta }{d\tau }=\frac{1}{8Q^2}\frac{\lambda }{2\pi a}_0^{2\pi }\frac{\chi a\mathrm{cos}^2\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}.$$
(33)
For $`p1`$ Eq. (33) can be written as,

$$\frac{d\vartheta }{d\tau }=\frac{1}{8Q^2}\frac{2\lambda }{\pi a}\left[1\frac{p^2}{4}\left(2\mathrm{ln}\frac{4}{p}1\right)\right].$$
(34)
The last term in the right-hand side of Eq. (34) describes a change of the frequency of small oscillations of the cantilever caused by the adiabatic reversals of $`\stackrel{}{M}`$. From Eq. (32) we have $`a(\tau )=a(0)\mathrm{exp}[\tau /(2Q)]`$. One can see from Eq. (34) that for the initial conditions (31) the influence of $`\stackrel{}{M}`$ results in decrease of the frequency of small oscillations of the cantilever in comparison with the case $`\lambda =0`$. For small $`p`$ (the value of $`p=\epsilon /(a\chi )\mathrm{exp}(t/2Q)`$ increases with time) the frequency of oscillations decreases when time increases, as shown in Fig. 5, while in the absence of the sample this frequency remains independent of time. We should note that in the studied approximation the sample does not influence the amplitude of the cantilever oscillations.

## VII Summary

We have studied theoretically and numerically the stationary cantilever vibrations in the novel OSCAR MRFM technique. Our results are based on the application of the classical theory for the motion of the cantilever and the paramagnetic moment of a cluster on the surface of the sample. We have estimated the resonant frequency shift for the cantilever vibrations. For the reasonable values of parameters our estimate is about 6 Hz per Bohr magneton. The sign of the shift depends on the initial direction of the paramagnetic moment relative to the initial position of the cantilever. We supported our estimation by the analytical analysis based on the perturbation theory and by the numerical solution of the equations of motion. Our perturbative approach is based on the fact that the influence of the paramagnetic moment on the sample is weak ($`\lambda 1`$). We considered also the regime of small oscillations of the paramagnetic moment near the transversal plane (linear OSCAR regime). Finally, we analyzed the damped oscillations of the cantilever (without the external force). We have shown that the frequency of the damped oscillations becomes time-dependent due to the adiabatic reversals of the paramagnetic moment.

## Acknowledgments

We are thankful to P.C. Hammel, D.V. Pelekhov, and D. Rugar for useful discussions. The work was supported by the Department of Energy (DOE) under contract W-7405-ENG-36, by the National Security Agency (NSA), by the Advanced Research and Development Activity (ARDA), and by the Defense Advanced Research Program Agency (DARPA) through the program MOSAIC.

## Appendix

Here we express the integral in Eq. (22) in terms of complete elliptic integrals,

$$_0^{\pi /2}\frac{\chi a\mathrm{cos}^2\psi d\psi }{\sqrt{\epsilon ^2+(\chi a\mathrm{cos}\psi )^2}}=_0^{\pi /2}\frac{\mathrm{cos}^2\psi d\psi }{\sqrt{p^2+\mathrm{cos}^2\psi }}=_0^{\pi /2}\frac{(1\mathrm{sin}^2\psi )d\psi }{\sqrt{p^2+1\mathrm{sin}^2\psi }}=$$

$$_0^{\pi /2}\frac{(p^2+1)\left(1\frac{1}{p^2+1}\mathrm{sin}^2\psi \right)p^2}{\sqrt{p^2+1}\sqrt{1\frac{1}{p^2+1}\mathrm{sin}^2\psi }}𝑑\psi ,$$

where we introduced the notation $`p=\epsilon /(a\chi )`$. Splitting this integral in two parts we obtain the right-hand side of Eq. (23).