text:
# Adiabatic processes need not correspond to optimal work

## Abstract

The minimum work principle states that work done on a thermally isolated equilibrium system is minimal for the adiabatically slow (reversible) realization of a given process. This principle, one of the formulations of the second law, is studied here for finite (possibly large) quantum systems interacting with macroscopic sources of work. It is shown to be valid as long as the adiabatic energy levels do not cross. If level crossing does occur, counter examples are discussed, showing that the minimum work principle can be violated and that optimal processes are neither adiabatically slow nor reversible.

The second law of thermodynamics landau ; balian ; perrot , formulated nearly one and half century ago, continues to be under scrutinity QL2L ; AN ; NA . While its status within equilibrium thermodynamics and statistical physics is by now well-settled landau ; balian ; perrot , its fate in various border situations is far from being clear. In the macroscopic realm the second law is a set of equivalent statements concerning quantities such as entropy, heat, work, etc. However, in more general situations these statements need not be equivalent and some, e.g. those involving entropy, may have only a limited applicabilityAN ; NA . In contrast to entropy, the concept of work has a well-defined operational meaning for finite systems interacting with macroscopic work sources balian . It is, perhaps, not accidental that Thomson’s formulation of the second law landau ; balian ; perrot  — no work can be extracted from an equilibrium system by means of a cyclis process — was proven ANThomson ; bassetttasaki  both in quantum and classical situation.

Here we study the minimum work principle which extends the Thomson’s formulation to non-cyclic processes landau ; balian ; perrot ; bassetttasaki , and provides a recipe for reducing energy costs. After formulating the principle and discussing it for macroscopic systems, we investigate it for finite systems coupled to macroscopic sources of work. Its domain of validity there is found to be large but definitely limited. These limits are illustrated via counterexamples.

The setup. Consider a quantum system S which is thermally isolated landau ; balian ; perrot : it moves according to its own dynamics and interacts with an external macroscopic work source. This interaction is realized via time-dependence of some parameters $`R(t)=\{R_1(t),R_2(t),\mathrm{}\}`$ of the system’s Hamiltonian $`H(t)=H\{R(t)\}`$. They move along a certain trajectory $`R(t)`$ which at some initial time $`t_\mathrm{i}`$ starts from $`R_\mathrm{i}=R(t_\mathrm{i})`$, and ends at $`R_\mathrm{f}=R(t_\mathrm{f})`$. The initial and final values of the Hamiltonian are $`H_\mathrm{i}=H\{R_\mathrm{i}\}`$ and $`H_\mathrm{f}=H\{R_\mathrm{f}\}`$, respectively. Initially S is assumed to be in equilibrium at temperature $`T=1/\beta 0`$, that is, S is described by a gibbsian density operator:

$`\rho (t_\mathrm{i})=\mathrm{exp}(\beta H_\mathrm{i})/Z_\mathrm{i},Z_\mathrm{i}=\mathrm{tr}e^{\beta H_\mathrm{i}}.`$ (1)
As usual, this equilibrium state is prepared by a weak interaction between S and a macroscopic thermal bath at temperature $`T`$ landau ; balian ; NA , and then decoupling S from the bath in order to achieve a thermally isolated process landau ; balian ; perrot .

The Hamiltonian $`H(t)`$ generates a unitary evolution:

$`i\mathrm{}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\rho (t)=[H(t),\rho (t)],\rho (t)=U(t)\rho (t_\mathrm{i})U^{}(t),`$ (2)
with time-ordered $`U(t)=\stackrel{}{\mathrm{exp}}[\frac{i}{\mathrm{}}_{t_\mathrm{i}}^tdsH(s)]`$. The work $`W`$ done on S reads landau ; balian ; perrot

$`W={\displaystyle _{t_\mathrm{i}}^{t_\mathrm{f}}}dt\mathrm{tr}[\rho (t)\dot{H}(t)]=\mathrm{tr}[H_\mathrm{f}\rho (t_\mathrm{f})]\mathrm{tr}[H_\mathrm{i}\rho (t_\mathrm{i})],`$ (3)
where we performed partial integration and inserted (2). This is the energy increase of S, which coincides with the energy decrease of the source.

The principle. Let S start in the state (1), and let $`R`$ move between $`R_\mathrm{i}`$ and $`R_\mathrm{f}`$ along a trajectory $`R(t)`$. The work done on S is $`W`$. Consider the adiabatically slow realization of this process: $`R`$ proceeds between the same values $`R_\mathrm{i}`$ and $`R_\mathrm{f}`$ and along the same trajectory, but now with a homogeneously vanishing velocity, thereby taking a very long time $`t_\mathrm{f}t_\mathrm{i}`$, at the cost of an amount work $`\stackrel{~}{W}`$. The minimum-work principle then asserts landau ; balian ; perrot

$`W\stackrel{~}{W}.`$ (4)
This is a statement on optimality: if work has to be extracted from S, $`W`$ is negative, and to make it as negative as possible one proceeds with very slow velocity. If during some operation work has to be added ($`W>0`$) to S, one wishes to minimize its amount, and operates slowly. For thermally isolated systems, adiabatically slow processes are reversible. This is standard if S is macroscopic landau ; balian ; perrot , and below it is shown to hold for a finite S as well, where the definition of reversibility extends unambiguously (i.e., without invoking entropy) perrot .

In macroscopic thermodynamics the minimum work principle is derived landau ; perrot  from certain axioms which ensure that, within the domain of their applicability, this principle is equivalent to other formulations of the second law. Derivations in the context of statistical thermodynamics are presented in jar ; narnhofer ; bassetttasaki . We discuss one of them now.

The minimal work principle for macroscopic systems is proven in two steps: first one considers the relative entropy $`\mathrm{tr}[\rho (t_\mathrm{f})\mathrm{ln}\rho (t_\mathrm{f})\rho (t_\mathrm{f})\mathrm{ln}\rho _{\mathrm{eq}}(H_\mathrm{f})]`$ between the final state $`\rho (t_\mathrm{f})`$ given by (2) and an equilibrium state $`\rho _{\mathrm{eq}}(H_\mathrm{f})=\mathrm{exp}(\beta H_\mathrm{f})/Z_\mathrm{f}`$, $`Z_\mathrm{f}=\mathrm{tr}e^{\beta H_\mathrm{f}}`$, a state corresponding to the final Hamiltonian $`H_\mathrm{f}`$ and the same temperature $`T=1/\beta `$. As follows from (2), $`\mathrm{tr}[\rho (t_\mathrm{f})\mathrm{ln}\rho (t_\mathrm{f})]=\mathrm{tr}[\rho (t_\mathrm{i})\mathrm{ln}\rho (t_\mathrm{i})]`$. This combined with (1, 3) and the non-negativity of relative entropy balian  yields:

$`WF(H_\mathrm{f})F(H_\mathrm{i})T\mathrm{ln}\mathrm{tr}e^{\beta H_\mathrm{i}}T\mathrm{ln}\mathrm{tr}e^{\beta H_\mathrm{f}},`$ (5)
where $`F(H_\mathrm{i})`$ and $`F(H_\mathrm{f})`$ are the gibbsian free energies corresponding to $`\rho (t_\mathrm{i})`$ and $`\rho _{\mathrm{eq}}(H_\mathrm{f})`$, respectively.

There are several classes of macroscopic systems for which one can show that the free energy difference in the RHS of (5) indeed coincides with the adiabatic work narnhofer ; jar ; NA .

Finite systems. For an arbitrary $`N`$-level quantum system S, Eq. (5) does not have the needed physical meaning, since in general $`F(H_\mathrm{f})F(H_\mathrm{i})`$ does not coincide with the the adiabatic work. It is knownflu  that for finite systems the final density matrix $`\rho (t_\mathrm{f})`$ given by (2) need not coincide with $`\rho _{\mathrm{eq}}(H_\mathrm{f})=\mathrm{exp}(\beta H_\mathrm{f})/Z_\mathrm{f}`$. This fact was recently applied for certain irreversible processesjapan .

Thus we need an independent derivation of (4). Let the spectral resolution of $`H(t)`$ and $`\rho (t_\mathrm{i})`$ be

$`H(t)={\displaystyle \underset{k=1}{\overset{N}{}}}\epsilon _k(t)|k,tk,t|,k,t|n,t=\delta _{kn},`$ (6)
$`\rho (t_\mathrm{i})={\displaystyle \underset{k=1}{\overset{N}{}}}p_k|k,t_\mathrm{i}k,t_\mathrm{i}|,p_k={\displaystyle \frac{e^{\beta \epsilon _k(t_\mathrm{i})}}{_ne^{\beta \epsilon _n(t_\mathrm{i})}}}.`$ (7)
At $`t=t_\mathrm{i}`$ we order the spectrum as

$`\epsilon _1(t_\mathrm{i})\mathrm{}\epsilon _N(t_\mathrm{i})p_1\mathrm{}p_N.`$ (8)
For $`t_\mathrm{i}tt_\mathrm{f}`$ we expand on the complete set $`|n,t`$:

$`U(t)|k,t_\mathrm{i}={\displaystyle \underset{n=1}{\overset{N}{}}}a_{kn}(t)e^{\frac{i}{\mathrm{}}_{t_\mathrm{i}}^tdt^{}\epsilon _n(t^{})}|n,t,`$ (9)
and use (3) to obtain:

$`W={\displaystyle \underset{k,n=1}{\overset{N}{}}}|a_{kn}(t_\mathrm{f})|^2p_k\epsilon _n(t_\mathrm{f}){\displaystyle \underset{k=1}{\overset{N}{}}}p_k\epsilon _k(t_\mathrm{i}).`$ (10)
A similar formula can be derived to express the adiabatic work $`\stackrel{~}{W}`$ in coefficients $`\stackrel{~}{a}_{kn}(t_\mathrm{f})`$. From the definition $`|a_{kn}(t_\mathrm{f})|^2=|n,t_\mathrm{f}|U|k,t_\mathrm{i}|^2`$ it follows that

$`{\displaystyle \underset{k=1}{\overset{N}{}}}|a_{kn}(t_\mathrm{f})|^2={\displaystyle \underset{n=1}{\overset{N}{}}}|a_{kn}(t_\mathrm{f})|^2=1.`$ (11)
With help of the identity: $`_{n=1}^N\epsilon _nx_n=\epsilon _N_{n=1}^Nx_n_{m=1}^{N1}[\epsilon _{m+1}\epsilon _m]_{n=1}^mx_n`$, we obtain using (10, 11) the general formula for the difference between non-adiabatic and adiabatic work:

$`W\stackrel{~}{W}=`$ $`{\displaystyle \underset{m=1}{\overset{N1}{}}}[\epsilon _{m+1}(t_\mathrm{f})\epsilon _m(t_\mathrm{f})]\mathrm{\Theta }_m,`$ (12)
$`\mathrm{\Theta }_m`$ $`{\displaystyle \underset{n=1}{\overset{m}{}}}{\displaystyle \underset{k=1}{\overset{N}{}}}p_k(|\stackrel{~}{a}_{kn}(t_\mathrm{f})|^2|a_{kn}(t_\mathrm{f})|^2).`$ (13)
To understand the meaning of this formula, let us first assume that the ordering (8) is kept at $`t=t_\mathrm{f}`$:

$`\epsilon _1(t_\mathrm{f})\mathrm{}\epsilon _N(t_\mathrm{f}).`$ (14)
If different energy levels did not cross each other (and equal ones do not become different), Eq. (14) is implied by Eq. (8). According to non-crossing rule mead1 , if $`H\{R\}`$ is real and only one of its parameters is varied with time, (14) is satisfied for any discrete-level quantum system: level-crossing, even if it happens in model-dependent calculations or due to approximate symmetry, does not survive arbitrary small perturbation where it is substituted by avoided crossing (for a more general $`H\{R\}`$ the conditions prohibiting level-crossing are more restrictive; see mead1 ). No level-crossings and natural conditions of smoothmess of $`H(t)`$ are sufficient for the standard quantum adiabatic theorem standardad  to ensure

$`\stackrel{~}{a}_{kn}(t_\mathrm{f})=\delta _{kn}.`$ (15)
Combined with (8, 15), Eq.(13) brings

$`\mathrm{\Theta }_m={\displaystyle \underset{k=1}{\overset{m}{}}}p_k(1{\displaystyle \underset{n=1}{\overset{m}{}}}|a_{kn}(t_\mathrm{f})|^2){\displaystyle \underset{n=1}{\overset{m}{}}}{\displaystyle \underset{k=m+1}{\overset{N}{}}}p_k|a_{kn}(t_\mathrm{f})|^2`$
$`p_m\left[m{\displaystyle \underset{k=1}{\overset{m}{}}}{\displaystyle \underset{n=1}{\overset{m}{}}}|a_{kn}(t_\mathrm{f})|^2{\displaystyle \underset{n=1}{\overset{m}{}}}{\displaystyle \underset{k=m+1}{\overset{N}{}}}|a_{kn}(t_\mathrm{f})|^2\right]=0.`$
Eqs. (10, 12, 14) together with $`\mathrm{\Theta }_m0`$ extend the minimum work principle (4) to cases where the adiabatic work is not equal to the difference in free energies.

Level crossing. The above non-crossing condition raises the question: Is the minimum work principle also valid if the adiabatic energy levels cross? Before addressing this question in detail, let us mention some popular misconceptions which surround the level-crossing problem: 1) The no-crossing rule is said to exclude all crossings. This is incorrect as the exclusion concerns situations where, in particular, only one independent parameter of a real Hamiltonian $`H\{R\}`$ is varied mead1 . Two parameters can produce robust level-crossing for such Hamiltonians. 2) It is believed that once levels can cross, $`\mathrm{\Delta }\epsilon 0`$, the very point of the adiabatic theorem disappears as the internal characteristic time $`\mathrm{}/\mathrm{\Delta }\epsilon `$ of S is infinite. This view misidentifies the proper internal time as seen below; see also avron  in this context. 3) It is sometimes believed that crossing is automatically followed by a population inversion. We shall find no support for that.

As a first example we consider a spin-$`1/2`$ particle with Hamiltonian

$`H(t)=h_1(s)\sigma _1h_3(s)\sigma _3,s=t/\tau ,`$ (16)
where $`\sigma _1`$, $`\sigma _3`$ and $`\sigma _2=i\sigma _1\sigma _3`$ are Pauli matrices, and where $`s`$ is the reduced time with $`\tau `$ being the characteristic time-scale. The magnetic fields $`h_1`$ and $`h_3`$ smoothly vary in time. Assume that i) for $`ss_\mathrm{i}<0`$ and for $`ss_\mathrm{f}>0`$, $`h_1(s)`$ and $`h_3(s)`$ go to constant values sufficiently fast; ii) for $`s0`$ one has: $`h_1(s)\alpha _1s^2`$, $`h_3(s)\alpha _3s`$, where $`\alpha _1`$ and $`\alpha _3`$ are positive constants. iii) $`h_1(s)`$ and $`h_3(s)`$ are non-zero for all $`s`$, $`s_\mathrm{i}ss_\mathrm{f}`$, except $`s=0`$. Not all these points are needed, but we choose them for clarity. Generalizations are indicated below. One writes (16) as

$`H=\left(\begin{array}{cc}\hfill h_3(s)& \hfill h_1(s)\\ \hfill h_1(s)& \hfill h_3(s)\end{array}\right)=\epsilon _1(s)\left(\begin{array}{cc}\hfill \mathrm{cos}\theta (s)& \hfill \mathrm{sin}\theta (s)\\ \hfill \mathrm{sin}\theta (s)& \hfill \mathrm{cos}\theta (s)\end{array}\right),`$ (20)
where $`\theta (s)=\mathrm{arctan}[h_1(s)/h_3(s)]`$ is a parameter in the interval $`\pi /2<\theta <\pi /2`$, and where

$`\epsilon _1(s)=\mathrm{sg}(s)\sqrt{h_3^2(s)+h_1^2(s)},\epsilon _2(s)=\epsilon _1(s)`$ (21)
are the adiabatic energy levels which cross at $`s=\theta (s)=0`$ ($`\sqrt{\mathrm{}}`$ is defined to be always positive). The adiabatic eigenvectors are, $`H(s)|k,s=\epsilon _k(s)|k,s`$, $`k=1,2`$,

$`|1,s=\left(\begin{array}{c}\hfill \mathrm{cos}\frac{1}{2}\theta (s)\\ \\ \hfill \mathrm{sin}\frac{1}{2}\theta (s)\end{array}\right),|2,s=\left(\begin{array}{c}\hfill \mathrm{sin}\frac{1}{2}\theta (s)\\ \\ \hfill \mathrm{cos}\frac{1}{2}\theta (s)\end{array}\right).`$ (28)
Both the eigenvalues and the eigenvectors are smooth functions of $`s`$. Eq. (14) is not valid, and (1013) imply:

$`W\stackrel{~}{W}=2\sqrt{h_1^2(s_\mathrm{f})+h_3^2(s_\mathrm{f})}\mathrm{\Theta }_1,\tau s_\mathrm{f}=t_\mathrm{f}.`$ (29)
Naively this already proves the violation. More carefully, our strategy is now to confirm (15) in the adiabatic limit $`\tau \mathrm{}`$ and thus to confirm that $`\mathrm{\Theta }_1>0`$, implying that the minimum work principle is indeed violated. To this end we apply the standard adiabatic perturbation theory standardad . Substituting (9) into (2) one has:

$`\dot{a}_{kn}={\displaystyle \underset{m=1}{\overset{N}{}}}a_{km}(t)e^{\frac{i}{\mathrm{}}_{t_\mathrm{i}}^tdt^{}[\epsilon _n(t^{})\epsilon _m(t^{})]}n|_t|m.`$ (30)
As $`|1`$ and $`|2`$ in (28) are real, $`n|n=1`$ implies $`n|_t|n=0`$. Since $`n|_t|m=\frac{1}{\tau }n|_s|m`$ the RHS of (30) contains a small parameter $`1/\tau `$. It is more convenient to introduce new variables: $`a_{kn}(t)=\delta _{kn}+b_{kn}(t)`$, $`b_{kn}(t_\mathrm{i})=0`$. To leading order in $`1/\tau `$, $`b_{kn}`$ can be neglected in the RHS of (30), yielding for $`a_{kn}(t_\mathrm{f})=b_{kn}(t_\mathrm{f})`$:

$`|a_{kn}(t_\mathrm{f})|^2=\left|{\displaystyle _{s_\mathrm{i}}^{s_\mathrm{f}}}dse^{\frac{i\tau }{\mathrm{}}_{s_\mathrm{i}}^sdu[\epsilon _k(u)\epsilon _n(u)]}n|_s|k\right|^2,`$ (31)
while $`|a_{kk}(t_\mathrm{f})|^2=1_{nk}|a_{kn}(t_\mathrm{f})|^2`$. In (31) we put $`s\tau =t`$, $`s^{}\tau =t^{}`$. For our model (1628), $`_{s_\mathrm{i}}^sdu[\epsilon _1(u)\epsilon _2(u)]=2_{s_\mathrm{i}}^sdu\epsilon _1(u)`$ has only one extremal point, at $`s=0`$. We also have from (28)

$`2|_s|1={\displaystyle \frac{\theta ^{}}{2}}={\displaystyle \frac{1}{2}}{\displaystyle \frac{h_1h_3^{}h_3h_1^{}}{h_3^2+h_1^2}},\theta ^{}{\displaystyle \frac{\mathrm{d}\theta }{\mathrm{d}s}}.`$ (32)
For large $`\tau `$ the integral in (31) can be calculated with use of the saddle-point method:

$`|a_{12}(t_\mathrm{f})|^2={\displaystyle \frac{\pi \mathrm{}}{\tau }}\left[{\displaystyle \frac{2|_s|1^2\sqrt{h_1^2+h_3^2}}{h_1h_1^{}+h_3h_3^{}}}\right]|_{s=0}={\displaystyle \frac{\pi \mathrm{}\alpha _1^2}{4\tau \alpha _3^3}}.`$ (33)
Eqs. (31, 33) extend the adiabatic theorem (15) for the level-crossing situation. More general statements of similar adiabatic theorems can be found in Ref. avron . Inserting $`\mathrm{\Theta }_1=(p_1p_2)|a_{12}(t_\mathrm{f})|^2>0`$ in Eq. (29) confirms the violation of the minimum work principle. Eq. (33) shows that the role of the proper internal characteristic time is played by $`\mathrm{}\alpha _1^2/\alpha _3^3`$ rather than by $`\mathrm{}/(\epsilon _1\epsilon _2)`$.

More generally, if $`\sqrt{h_3^2(s)+h_1^2(s)}`$ is a smooth function for all real $`s`$ (e.g., it is not $`|s|`$), there are no crossings of eigenvalues and (4) is valid. If both $`h_1(s)`$ and $`h_3(s)`$ are linear for $`s0`$, the leading term presented in (33) vanishes due to $`2|_s|1^2|_{s=0}=0`$, and one needs the second-order in the saddle-point expansion, to be compared with the second-order term of the adiabatic perturbation theory. This leads to the same physical conclusions as (33) did, but with $`|a_{12}(t_\mathrm{f})|^2\tau ^3`$.

One can calculate $`|a_{kn}(t_\mathrm{f})|`$ yet in another limiting case, where the characteristic time $`\tau `$ is very short. It is well-known standardad  that in this limit energy changes can be calculated with frozen initial state of S. For the present situation this leads from (28) to $`|a_{12}(t_\mathrm{f})|^2=|a_{21}(t_\mathrm{f})|^2=|1,t_\mathrm{f}|\mathrm{\hspace{0.17em}2},t_\mathrm{i}|^2=\mathrm{sin}^2\frac{1}{2}[\theta (t_\mathrm{f})\theta (t_\mathrm{i})]`$, and thus to $`\mathrm{\Theta }_1=(p_1p_2)\mathrm{sin}^2\frac{1}{2}[\theta (t_\mathrm{f})\theta (t_\mathrm{i})]`$, again positive.

Exactly solvable model with level crossing. Consider a two-level system with Hamiltonian

$`H(t)=\mathrm{}\omega \left(\begin{array}{cc}\hfill s\mathrm{cos}^2s& \hfill \frac{1}{2}s\mathrm{sin}2s\\ & \\ \hfill \frac{1}{2}s\mathrm{sin}2s& \hfill s\mathrm{sin}^2s\end{array}\right),s={\displaystyle \frac{t}{\tau }},`$ (37)
where $`\tau `$ is the characteristic time-scale, and $`\omega `$ is a constant. For $`s_\mathrm{i}<0`$ denote the adiabatic energy levels as $`\epsilon _1(s_\mathrm{i})=\mathrm{}\omega s_\mathrm{i}<\epsilon _2(s_\mathrm{i})=0`$. They cross when $`s`$ passes through zero. Eq. (30) for this model can be solved exactly in terms of hypergeometric functions. Postponing the detailed discussion, we present in Fig. 1 the behavior of $`|a_{12}(s_\mathrm{f})|^2`$ as a function of $`\tau `$.

Since from (12, 37) one has for the present model $`W\stackrel{~}{W}=\mathrm{}\omega s_\mathrm{f}|a_{12}(s_\mathrm{f})|^2\mathrm{tanh}(\frac{1}{2}\beta \mathrm{}\omega |s_\mathrm{i}|)`$, violations of the minimum work principle exist for $`s_\mathrm{f}>0`$, and they are maximal for $`|a_{12}(s_\mathrm{f})|^21`$. This is seen to be the case in Fig. 1 for some $`\tau `$ near $`\tau =1`$. Notice that both the first-order perturbation theory and the saddle-point approximation are adequately reproducing $`|a_{12}(s_\mathrm{f})|^2`$ for $`\tau 10`$.

Let S has a finite amount of levels, and two of them cross. For quasi-adiabatic processes ($`\tau `$ is large but finite) the transition probability between non-crossing levels is exponentially small standardad ; pechu , while as we saw it has power-law smallness for two crossing levels. Then one neglects in (10) the factors $`|a_{kn}(t_\mathrm{f})|^2`$ coming from any non-crossed levels $`k`$ and $`n`$, and the problem is reduced to the two-level situation. Thus already one crossing suffices to detect limits of the minimum work principle. The reduction to the two-level situation takes place also in a macroscopic system which has few discrete levels at the bottom of a continuous spectrum, since for low temperatures these levels can decouple from the rest of the spectrum.

Cyclic processes and reversibility. The above results do not imply any violation of the second law in Thomson’s formulationANThomson : no work is extracted from S during a cyclic process, $`W_\mathrm{c}0`$. We illustrate its general proof in the context of the level crossing model given by (1628). Assume that the trajectory $`R(t)=(h_1(t),h_2(t))`$ described there is supplemented by another trajectory $`R^{}(t)`$ which brings the parameters back to their initial values $`(h_1(t_\mathrm{i}),h_3(t_\mathrm{i}))`$ so that the overall process $`R+R^{}`$ is cyclic. If $`R^{}`$ crosses the levels backwards, then at the final time of $`R^{}`$ Eq. (14) is valid, and (1015) imply:

$`W_\mathrm{c}=|a_{12}|^2(p_1p_2)[\epsilon _2(t_\mathrm{i})\epsilon _1(t_\mathrm{i})]\stackrel{~}{W}_\mathrm{c}=0,`$ (38)
where $`|a_{12}|^21`$ now corresponds to the full cyclic process $`R+R^{}`$. Eq. (38) confirms the intuitive expectation that non-adiabatic process are less optimal. In particular, this is valid if $`R^{}`$ is exactly the same process $`R`$ moved backwards with the same speed. Then $`\stackrel{~}{W}_\mathrm{c}=0`$ means that $`R`$ is a reversible process in the standard thermodynamical sense landau ; balian ; perrot . If $`R^{}`$ does not induce another level crossing, i.e., $`h_1(s)`$ and $`h_2(s)`$ in Eq. (16) return to their initial values without simultaneously crossing zero, then $`\epsilon _1(t_\mathrm{f})=\epsilon _2(t_\mathrm{i})`$, $`\epsilon _2(t_\mathrm{f})=\epsilon _1(t_\mathrm{i})`$ and Eqs. (10, 15) imply

$`\stackrel{~}{W}_\mathrm{c}=(p_1p_2)[\epsilon _2(t_\mathrm{i})\epsilon _1(t_\mathrm{i})]W_\mathrm{c}=|a_{11}|^2\stackrel{~}{W}_\mathrm{c}>0.`$
In contrast to (38), non-adiabatic processes are more optimal if $`R+R^{}`$ contains one level-crossing (or an odd number of them). We thus have found here a violation of the minimum work principle for a cyclic process.
In conclusion, we have studied the minimum work principle for finite systems coupled to external sources of work. As compared to other formulations of the second law, this principle has a direct practical meaning as it provides a recipe for reducing energy costs of various processes. We gave its general proof and have shown that it may become limited if there are crossings of adiabatic energy levels: optimal processes need to be neither slow nor reversible. Already one crossing suffices to note violations of the principle. If this is the case, the optimal process occurs at some finite speed.
Level-crossing was observed, e.g., in molecular and chemical physics chem . It is not a rare effectrare : if the number of externally varied parameters is larger then two, then for typical spectra level crossings are even much more frequent than avoided crossingsrare . It is possible that the presented limits of the minimum work principe may serve as a test for level crossings.
Together with the universal validity of Thomson’s formulation of the second law ANThomson ; bassetttasaki , the limits of the principle imply that the very equivalence between various formulations of the second law may be broken for a finite system coupled to macroscopic sources of work: different formulations are based on different physical mechanisms and have different ranges of validity. Similar results on non-equivalence of various formulations of the second law were found in Ref. AN ; NA , where for a quantum particle coupled to a macroscopic thermal bath, it was shown that some formulations, e.g., Clausius inequality and positivity of the energy dispersion rate, are satisfied at sufficiently high temperatures of the bath, but can be invalid at low temperatures, that is, in the quantum regime.
The work is supported by FOM/NWO (Netherlands).


format: text

text:
# Moving and staying together without a leader

## 1 Introduction

The emergence of collective motion of self-propelled organisms (bird flocks, fish schools, herds, slime molds, bacteria colonies, etc.) is a fascinating phenomenon which attracted the attention of (theoretical) physicists only recently (). Particularly intriguing are the situations where no “leader” with specific properties is present in the group, no mediating field helps organizing the collective dynamics (e.g. no chemotaxis), and interactions are short-range. In this case, even the possibility of collective motion may seem surprising.

However “simple” the involved organisms may be, they are still tremendously complex for a physicist and his inclination will often be to go away from the detailed, intricate, as-faithful-as-possible modeling approach usually taken by biologists, and to adopt “minimal models” hopefully catching the crucial, universal properties which may underlie seemingly different situations.

In this setting, the organisms can be reduced to points which move at finite velocity and interact with neighbors. This is in fact what Vicsek and collaborators did when introducing their minimal model for collective motion.

## 2 Vicsek’s model

Vicsek’s model  consists in pointwise particles labeled by $`i`$ which move synchronously at discrete timesteps $`\mathrm{\Delta }t`$ by a fixed distance $`v_0`$ along a direction $`\theta _i`$. This angle is calculated from the current velocities of all particles $`j`$ within an interaction range $`r_0`$, reflecting the only “force” at play, a tendency to align with neighboring particles:

$$\theta _i^{t+1}=\mathrm{arg}\left[\underset{ji}{}\stackrel{}{v}_j^t\right]+\eta \xi _i^t,$$
(1)
where $`\stackrel{}{v}_i^t`$ is the velocity vector of magnitude $`v_0`$ along direction $`\theta _i`$ and $`\xi _i^t`$ is a delta-correlated white noise ($`\xi [\pi ,\pi ]`$). Fixing $`r_0=1`$, $`\mathrm{\Delta }t=1`$, and choosing, without loss of generality, a value $`v_0\mathrm{\Delta }t<r_0`$, Vicsek et al studied the behavior of this simple model in the two dimensional parameter space formed by the noise strength $`\eta `$ and $`\rho `$, the particle density. They found, at large $`\rho `$ and/or small $`\eta `$, the existence of an ordered phase characterized by:

$$V|\stackrel{}{v}_i^t_i|_t>0,$$
(2)
i.e. a domain of parameter space in which the particles move collectively.

The existence of the ordered phase was later proved analytically  via a continuous model for the coarse-grained particle velocity and density. Vicsek et al devoted most of their effort to studying the transition to the ordered phase . They found numerically a continuous transition characterized by scaling laws and they tried to estimate the corresponding set of critical exponents.

## 3 Collective and cohesive motion

Vicsek’s model accounts rather well, at least at a qualitative level, for situations where the organisms interact at short distances but need not stay together. This is for instance the case of the bacterial bath recently studied by Wu and Libchaber .

In this experiment, E. Coli bacteria are swimming freely within a fluid film of thickness approximately equal to their size. By seeding the system with polystyrene beads and recording the trajectories of these passive tracers, Wu and Libchaber showed that the bacteria perform superdiffusive motion crossing over to normal diffusion. We later argued that the superdiffusive behavior is likely to be due to the onset of collective motion as in Vicsek’s model .

When the situation to be described involves the overall cohesion of the population, Vicsek’s model needs to be supplemented by a suitable feature. Indeed, an initially cohesive flock of particles will disperse in an open space. In other words, no collective motion is possible in the zero density limit of this model. In the following, we extend Vicsek’s model to account for the possible cohesion of the population of particles.

The above remark is by no means new. Early models for the collective motion of “boids” (contraction of “birdoid”, a term used by computer animation graphics specialists) do include a two-body repulsive-attractive interaction . More recent works by physicists also included this ingredient, but they either comprised an extra global interaction , or the actual interaction range used extended over the whole flock for the sizes considered, making it effectively global . Another encountered pitfall, from our point of view at least, is to enforce the cohesion by the confinement to a rather small, close, space .

Here we want to be, in a sense, in the least-favorable circumstances for observing collective motion: no leader in the group, strongly noisy environment and/or communications, strictly local interactions, and no confinement at all. The “minimal” model presented below is one of the simplest possible ones satisfying these constraints.

## 4 A minimal model

In addition to the possibility of achieving cohesion, we also want to confer a “physical” extent to the particles, a feature absent from Vicsek’s point-particles approach. Adding a Lennard-Jones-type body force $`\stackrel{}{f}`$ acting between each pair of particles within distance $`r_0`$ from each other offers such a possibility.

Equation (1) is then replaced by

$$\theta _i^{t+1}=\mathrm{arg}\left[\alpha \underset{ji}{}\stackrel{}{v}_j^t+\beta \underset{ji}{}\stackrel{}{f}_{ij}\right]+\eta \xi _i^t,$$
(3)
where $`\alpha `$ and $`\beta `$ control the relative importance of the two “forces”. The precise form of the dependence of the body force on the distance between the two particles involved is not important. It is enough to ensure a hard-core repulsion at distance $`r_\mathrm{c}`$ and an “equilibrium” preferred distance $`r_\mathrm{e}`$.

In the following, we use

$$\stackrel{}{f}_{ij}=\stackrel{}{e}_{ij}\{\begin{array}{cc}\mathrm{}\hfill & \mathrm{if}r_{ij}<r_\mathrm{c},\hfill \\ \frac{1}{4}\frac{r_{ij}r_\mathrm{e}}{r_\mathrm{a}r_\mathrm{e}}\hfill & \mathrm{if}r_\mathrm{c}<r_{ij}<r_\mathrm{a},\hfill \\ 1\hfill & \mathrm{if}r_\mathrm{a}<r_{ij}<r_0.\hfill \end{array}$$
(4)
with $`r_{ij}`$ the distance between boids $`i`$ and $`j`$, $`\stackrel{}{e}_{ij}`$ the unit vector along the segment going from $`i`$ to $`j`$, and the numerical values $`r_\mathrm{c}=0.2`$, $`r_\mathrm{e}=0.5`$ and $`r_\mathrm{a}=0.8`$.

We have also tested other types of noise term in the model. In particular, considering the noise as the uncertainty with which each boid “evaluates” the force exerted on itself by the $`N_i`$ neighboring boids leads to change Eq. (3) to:

$$\theta _i^{t+1}=\mathrm{arg}\left[\alpha \underset{ji}{}\stackrel{}{v}_j^t+\beta \underset{ji}{}\stackrel{}{f}_{ij}+N_i\eta \stackrel{}{u}_i^t\right]$$
(5)
where $`\stackrel{}{u}_i^t`$ is a unit vector of random orientation. There are delicate issues related to the choice of the noise term, in particular with respect to the critical properties of the transition to collective motion . We mostly considered, in the following, the noise term as prescripted above in Eq. (5). A detailed analysis of the influence of the nature of the noise is left for future studies, but we are confident that the results presented here hold generally, at least at a qualitative level.

Finally, in order to ensure that each particle only interacts with its “first layer” of neighbors (within distance $`r_0`$), we calculate, at each timestep, the Voronoi tessellation of the population —this has the additional advantage of providing a natural definition of the “cells” associated with each particle. The interacting neighbors are then retricted to be those of the particles within distance $`r_0`$ which are also neighbors in the Voronoi sense.

## 5 Typical phases

One can easily guess the “phases” that the above model can exhibit for a fixed noise strength $`\eta `$ ($`\eta =1.0`$ in the following, for a summary of parameters see Table 1). We now present them in a qualitative manner. The results presented below were all obtained in the two-dimensional case, but most of them hold in three dimensions.

When the body force is weak (small $`\beta `$ values), the cohesion of a flock cannot be maintained. In a finite box (finite particle density $`\rho `$), one is left with a gas-like phase (not shown). An arbitrarily large flock in an infinite space disintegrates, eventually leaving isolated random-walking particles.

For large enough $`\beta `$, we can expect the cohesion to be maintained. Figure 1 shows such cohesive flocks. The internal structure of the flocks depends also on $`\beta `$: for large body force, positional quasi-order is present, and the particles locally form an hexagonal crystal (Fig. 1ab). For intermediate values of $`\beta `$, no positional order arises, and the flock behaves like a liquid droplet (Fig. 1cd).

The influence of the alignment “force” is manifested by the global motion of the flock: for large enough $`\alpha `$ values, the flock moves ($`V>0`$). Depending on $`\beta `$, one has then either a “moving droplet” (Fig. 1d) or a “flying crystal” (Fig. 1b).

## 6 Order parameters

Order parameters have to be defined to allow for a quantitative distinction between the phases described above.

The limit of cohesion separating the “liquid” phases from the “gas” can be determined by measuring the distribution of the sizes of particles clusters, thanks to an implementation of the Hoshen-Kopelman  algorithm. A cohesive flock is then one for which $`n`$, the size of the largest cluster, is of order $`N`$, the total number of particles. Below, we use the criterion $`n/N=\frac{1}{2}`$ to define the transition. Increasing $`\beta `$, $`n/N`$ sharply rises to order-one values (Fig. 2a).

The “liquid/solid” transition takes place when $`\beta `$ is large enough so that cohesion of the population is ensured. To determine this onset of positional quasi-order within a finite but arbitrarily large cohesive flock, we first observe that, whether in the collective motion region or not, the “liquid” and “solid” phases can be distinguished by the fact that particles diffuse with respect to each other in the “liquid”, whereas neighboring particles always remain close to each other in the “solid” (Fig. 3). To be more precise, in the “liquid”, initially close-by particles remain so for some trapping time $`\tau `$ (which can be defined or used to collapse the curves of Fig. 3). Approaching the “solid” phase (by increasing $`\beta `$), $`\tau `$ diverges. In addition, since we are not dealing with a translation-invariant system, we should also distinguish between the diffusion properties of the particles depending on their relative position within the flock. We have thus defined different “sectors” (“core”, “head”, “tail”, “sides”) as explained in Fig. 4.

If all sectors are roughly equivalent in non-moving droplets (at least sufficiently far from the “liquid/solid” transition), some differences are observed within moving droplets: the outer regions and in particular the head are more active, whereas cohesion is stronger in the core (Fig. 5). Thus, in principle, different trapping times $`\tau `$ can be defined for the different regions of a (moving) flock. But the depth of the outer, “more liquid”, layer does not depend on the flock size (if it is big enough), so that, in large flocks, most of the population behaves as the core. Consequently, the relative diffusion averaged over the whole flock suffers from finite-size effects, but they disappear in the large-size limit (see below). To sum up, the trapping time $`\tau `$ (measured on all particles of the flock) is a good quantity to track the “liquid/solid” transition. However, instead of directly estimating of $`\tau `$, we measured $`\mathrm{\Delta }`$, the relative diffusion over some large time $`T`$ of initially-neighboring particles:

$$\mathrm{\Delta }\frac{1}{n_i}\underset{ji}{}\left(1\frac{r_{ij}^2(t)}{r_{ij}^2(t+T)}\right)_{i,t},$$
(6)
where the $`n_i`$ particles $`j`$ are the neighbors of particle $`i`$ at time $`t`$. Time $`T`$ is taken to be proportional to the volume of the system, which ensures that $`\mathrm{\Delta }`$ records, in the large-size limit, an asymptotic property of the system (because $`T\tau `$). Clearly, $`\mathrm{\Delta }1`$ in the “liquid” phase for a large enough system, while $`\mathrm{\Delta }0`$ in the “solid” phase. Indeed, increasing $`\beta `$, $`\mathrm{\Delta }`$ falls off sharply (Fig. 2b). The transition point is chosen to be at $`\mathrm{\Delta }=\frac{1}{2}`$.

Finally, we must be able to distinguish the regimes for which collective motion arises. To this aim, we use the average velocity $`V`$, as defined in (2). Increasing $`\alpha `$, $`V`$ reaches order-one values. The onset of collective motion is chosen to be at $`V=v_0/2`$. In the gas phase, one expects $`V=0`$ independently of the strengh of the alignment force. Nevertheless, this phase is rather sharply divided in two when studying the model at finite particle density. The largest cluster size may then be small ($`nN`$), but it is (almost) always larger than one. At any given time, thus, $`n>1`$, and the collective motion order parameter $`V`$ can be defined restricted to the particles belonging to the largest cluster. (Since clusters merge and break, the particles involved generally change along time.)

## 7 Phase diagram

After a brief discussion of the nature of the transitions involved, we first present the phase diagram of our model for a finite density of particles in a large and fixed box size. Then we estimate finite-size effects on the location of the phase boundaries. Finally, we argue that the phase diagram can also be defined in the zero-density limit where an arbitrarily large flock wanders in an infinite space.

### 7.1 Nature of the transitions

As expected in usual phase transitions, we found that in our model the “gas/liquid” and “liquid/solid” transitions are first-order. In insets of Fig. 2a and b, we show the evolution of the probability distribution function of the order parameter as one crosses these transition lines. The bimodal character of these pdf at the transition is typical of first-order phase transitions, indicating the coexistence of two metastable states. At the “gas/liquid” transition point, dispersed and aggregated boids coexist and there are exchanges between the two phases along time. At the “liquid/solid” transition point, cohesion is ensured and one observes the quasi-frozen regions in an otherwise more “liquid” flock. The quasi-solid parts are often located in the core of the flock. (Fig. 6).

The nature of the transition for the onset of collective motion is a delicate issue in our model. Whereas its second-order character is rather well-established for Vicsek’s core model, we recently discovered that, in fact, the implementation of the noise term may change the nature of the transition. The detailed investigation of this, in particular in presence of the cohesive force, will be presented in a future publication .

### 7.2 Fixed population and fixed box size

A systematic scan of the ($`\alpha `$, $`\beta `$) parameter plane was performed for a flock of $`N=2025`$ particles living on a square surface of linear size $`L=180`$ ($`\rho =1/16`$) with periodic boundary conditions. Using the criteria defined above, we obtained the phase diagram presented in Fig. 7.

For each parameter value, we used an initially-aggregated flock. We let the system evolve during a time $`\tau L^d`$, and then we recorded each order parameter and its histogram along time. The transition points were determined by dichotomy, the precision of which is reflected in the error bars.

The basic expected features are found: the horizontal “gas/liquid” and “liquid/solid” transitions are crossed by the vertical “moving/non-moving” line. Near this line, however, one observes a strong deformation of the “gas/liquid” and “liquid/solid” boundaries. This cannot be understood without a careful study of the collective motion transition .

Note also that the “gas” phase itself is crossed by the line marking the onset of collective motion, using, as explained above, the average velocity $`V`$ of the $`n`$ particles of the largest cluster as the order parameter.

### 7.3 Finite size and saturated vapour effects

We are ultimately interested in the possibility of collective and cohesive motion for an arbitrarily-large flock in an infinite space. The phase diagram of Fig. 7 was obtained at a fixed system size and constant density. Thus both limits of infinite-size and zero-density have to be taken to reach the asymptotic regime of interest. Of course this is mostly relevant to the onset of cohesion (the “gas/liquid” transition). Here we first study each limit separately, i.e. we investigate finite-size effects at fixed particle density and expansion at fixed particle number. Then we discuss the double-limit regime of interest.

Performing such a task for the whole parameter plane far exceeds our available computer power. We restricted ourselves to three typical cases: in the non-moving phase ($`\alpha =1.0`$), in the moving phase ($`\alpha =3.0`$) and near the transition to collective motion ($`\alpha =1.75`$ for the “gas/liquid” transition and $`\alpha =2.1`$ for the “liquid/solid” transition).

In simulations performed at $`\rho =1/16`$, varying $`N`$ and $`L`$, the transition points $`\beta _{\mathrm{GL}}^\rho (N)`$ and $`\beta _{\mathrm{LS}}^\rho (N)`$ converge exponentially as $`L=\sqrt{N/\rho }`$ increases. In Fig. 8, we show this for both the “gas/liquid” (a) and the “liquid/solid” (b) transitions for $`\alpha =3`$ (moving phases). This allows to determine asymptotic transition points. Note that exponentially-decreasing finite-size effects are typical of first order phase transition at equilibrium . In the two other cases (non-moving phase, and onset of motion), the results are similar, with the asymptotic values being quickly reached in the non-moving phase ($`\alpha =1`$) .

In simulations performed at fixed $`N`$ varying $`L`$, we study instead the expansion of the system. The (finite) spatial extent induces a confinement effect which increases the pressure at the coexistence point. We thus expect a displacement of the transition point. We find that $`\beta _{\mathrm{GL}}(\rho )`$ also converges exponentially, and thus transition points $`\beta _{\mathrm{GL}}(N)`$ are well-defined. In Fig. 8c, we show this only for the “gas/liquid” transition at the $`\alpha `$ values used in Fig. 8a, since the “liquid/solid” transitions points were observed to be independent of the box size.

### 7.4 Zero-density limit

The above results provide evidence that our model possesses well-defined, asymptotic phase diagrams at either fixed particle density or at fixed number of particles. The double limit mentioned above can be approached in essentially three different ways.

One can take one limit after the other one, repeating either the calculations of Fig. 8a at lower and lower densities, or those of Fig. 8d at larger and larger flock size. However, this straightforward program involves very heavy numerical simulations. Therefore, we only considered two cases ($`\alpha =1.75`$ and $`\alpha =3.0`$). As expected, the transition points converge (exponentially) independently of the order with which the two limits are taken, yielding estimates of the zero-density limit (Fig. 9). These estimates are compatible with each other. There are three different sources of error: statistical and systematic errors when evaluating the order parameters of each system, and then fitting errors when determinating $`\beta _{\mathrm{GL}}(N)`$ and $`\beta _{\mathrm{GL}}(\rho )`$. For $`\alpha =1.75`$, we find a $`9\%`$ difference in the estimated asymptotic values whose origin is probably the statistical errors on the larger flocks simulations.

### 7.5 Evaporation of a flock

There exists a third manner of approaching the zero-density limit of the “gas/liquid” transition. It consists in quenching a cohesive flock observed at large-enough $`\beta `$ to a lower $`\beta `$ value. If the flock is quenched below the “gas/liquid” line, it will “evaporate”.

The largest cluster will progressively lose particles before finally equilibrating in the gas phase. This transient can be expected to be governed by an effective surface tension $`\gamma `$ and the boundary of the largest cluster should then be governed by some Allen-Cahn law : $`v_n=\gamma \kappa `$ where $`v_n`$ is the (local) normal velocity and $`\kappa `$ the local curvature. Assuming that the mass and the surface of the main cluster remain proportional, in the two-dimensional problem: $`n=\rho _{\mathrm{loc}}\pi R^2`$, we can integrate this equation, and we obtain the relation :

$$n=N\lambda t\mathrm{with}\lambda =\pi \rho _{\mathrm{loc}}\gamma $$

where $`\rho _{\mathrm{loc}}`$ is the local density: the size of a circular flock should decrease linearly in time, with the proportionality constant providing an estimate of $`\lambda `$.

This is indeed what can be observed in our model (Fig. 10a and b). In these experiments, an initially large cohesive flock is prepared at some $`\beta `$ value above the transition. The system size is taken so as be close to the zero-density limit ($`\rho =1/256`$ in Fig. 10). We measured $`T_{\mathrm{ev}}`$, the time taken by the largest cluster to reach a given normalized mass $`n/N`$, as well as its surface at the same mass/time. We thus checked the proportionality between mass and surface and between mass and time, after transients and before the system approaches the equilibrium (Fig. 10a and b). Note that due to the abrupt change of parameters the flock actually first expands after the quench before setting in the “true” evaporation regime (Fig. 10b). This experiment can, at first sight, be thought of being free of confinement effects, and the Allen-Cahn law is expected to be satisfied at all times.

Figure 10c shows the dependence of $`\lambda `$ on $`\beta `$: $`\lambda `$ quickly decays and reaches very small values for $`\beta 0.8`$. On a logarithmic scale (Fig. 10d), we can distinguish two exponential regimes, on each side of this value, which corresponds roughly to the finite-size threshold $`\beta _{\mathrm{GL}}(L)`$ determined above (and is rather far from the asymptotic threshold determined above to be around 1.5) A simple argument can account for this behavior: Consider a flock of $`N`$ boids and suppose there is no short-time expansion (such as seen in Fig. 10b). The mass of the largest cluster $`n(t)`$, from then on, decreases linearly with time until it reaches, at time $`T_{\mathrm{}}`$, the equilibrium value $`n_{\mathrm{eq}}`$. At every time $`t`$, we have:

$$Nn(t)=\frac{t}{T_{\mathrm{}}}(Nn_{\mathrm{eq}})=\lambda t.$$
(7)
From the theory of first order phase transitions of systems at equilibrium (see for instance ), we expect the order parameter $`n_{\mathrm{eq}}/N`$ to behave like

$$\frac{n_{\mathrm{eq}}}{N}\frac{1}{2}\left[1+\mathrm{tanh}\left(K_1L^d(\beta \beta _{\mathrm{GL}})\right)\right],$$

in the vicinity of the transition, where $`K_1`$ is a constant. Moreover, we expect that $`T_{\mathrm{}}`$ depends linearly on $`N`$ (see Fig. 10a, where results for two different system sizes have been superimposed).

From a mean-field point of view, $`T_{\mathrm{}}/N`$ can be interpreted as the mean time required for a particle to escape from the interaction of another boid. Given the interaction we use (see Eq. (4)), we can assume that the potential is harmonic, so that the escape time is proportional to the exponential of the potential depth $`\beta `$. Finally, we get

$`\lambda `$ $`=`$ $`{\displaystyle \frac{N}{T_{\mathrm{}}}}\left(1{\displaystyle \frac{n_{\mathrm{eq}}}{N}}\right),\mathrm{or}`$
$`\lambda `$ $``$ $`\mathrm{exp}\left(K_2\beta \right)\left[1\mathrm{tanh}\left[K_1L^d\left(\beta \beta _{\mathrm{GL}}(L)\right)\right]\right]`$
Approximating $`\mathrm{tanh}`$ by $`(1+\mathrm{exp})`$ and $`(1\mathrm{exp})`$ sufficiently far below and above $`\beta =\beta _{\mathrm{GL}}(L)`$, we find the two exponential regimes mentioned above. More precisely, for $`\beta <\beta _{\mathrm{GL}}(L)`$, the surface tension should be almost independent on system size, whereas for $`\beta >\beta _{\mathrm{GL}}(L)`$, $`\lambda `$ is governed by the finite size effects and its slope increases like $`L^d`$. Quantitative agreements with the above approximation would require too large numerical calculations, but our partial data is consistent with the above predictions.

A concluding remark to this investigation of the effective surface tension governing the evaporation of a flock is that, contrary to naive arguments, this method does not offer much advantage over the double-limit procedure presented in Section 7.4. Indeed, as shown above, it only allows a rather easy determination of $`\beta _{\mathrm{GL}}(L)`$, while the finite-size effects remain hard to estimate quantitatively.

## 8 Micro vs macro motion

The existence of cohesive phases being now well-established, a natural question is that of the properties of the trajectories of cohesive flocks (the “macroscopic” motion) and it is interesting to compare those to the trajectories of the individuals composing the flock (“microscopic” motion). Postponing again the account of what happens in this respect near the onset of collective motion to a further publication , we studied, for the four possible cohesive phases, the mean square displacement $`\mathrm{\Delta }r^2(t)`$ of the center of mass, as well as of individual boids.

Our model being essentially stochastic, it is no surprise that, at large times, we observe that $`\mathrm{\Delta }r^2(t)t`$, i.e. the flock performs Brownian motion, in all cases (Fig. 11 and 12). When in a moving phase (either “liquid” or “solid”), this random walk may consist of ballistic flights separated by less coherent intervals during which the flock often changes direction. This is testified by the ballistic part ($`\mathrm{\Delta }r^2(t)t^2`$) of the plots in Fig. 11b and by the trajectories itselves 12.

The moving and non-moving phases can also be distinguished by opposite finite-size effects: in the non-moving phases, the diffusion constant of the macroscopic random walk decreases with system size (like $`1/N`$), whereas the ballistic flights’ duration increases with system size in the moving phases (Fig. 13).

Comparing the macroscopic motion of “liquid” and “solid” flocks is not well defined since, in the phase diagram, the line marking the onset of collective motion is not straight. Nevertheless, at comparable distances from this line, one notices that ballistic flights tend to be longer for flying crystals than for moving droplets. Similarly, the diffusion constant of “solid”, non-moving flocks is smaller than that of non-moving droplets.

Finally, even though we have already considered the mutual dispersion of initially-neighboring boids (see Fig. 3), we also studied the diffusive properties of individual boids within cohesive flocks, substracting out the translation motion of the center of mass. In “solid” flocks, one can hardly record any such motion. For droplets (Fig. 14), on the other hand, this microscopic motion depends on the macroscopic motion. When the droplet is fixed, the diffusion is normal $`\mathrm{\Delta }r^2t`$, whereas a boid which belongs to moving flock diffuses as $`\mathrm{\Delta }r^2t^\alpha `$, with $`\alpha 4/3`$. We interpret this as being due to some mesoscopic “hydrodynamical” structures within moving flocks (jets, vortices, etc.). J. Toner and Y. Tu  have shown via a mesoscopic equation that collective motion induces transverse correlations even in a co-moving frame. Therefore, boid diffusion must be faster than Brownian. They predicted an exponent equal to $`4/3`$, with which our results are in good agreement (Fig. 14, top line).

## 9 Summary and Perspectives

We have introduced a simple model for the collective and cohesive motion of self-propelled particles. We have described its various dynamical phases, defined order parameters to distinguish them, and presented a typical phase diagram at large but finite number of particles $`N`$ and large but finite system size $`L`$ (Fig. 7). Even though we have provided evidence that this phase diagram possesses well-defined $`N\mathrm{}`$ and $`L\mathrm{}`$ limits, these limit diagrams require too heavy numerical simulations to be determined at this stage.

We have also argued that the double limit of an arbitrarily large flock evolving in an infinite space is also well-defined. Here also, the mapping of the phase diagram in this limit is currently out of reach. But the existence of cohesive phases in a model of noisy, short-range interaction, identical particles is ensured, which was one of our primary goals. Futhermore, we can sketch the zero-density asymptotic diagram from our partial knowledge (Fig. 15).

A few remarks are in order to explain the expected shape of the asymptotic “gas/liquid” boundary: in the non-moving phase, this transition is almost independent on $`\alpha `$ and size-effects are negligible. We thus expect this line to be horizontal, in agreement with mean-field arguments (see Appendix). Similarly, our preliminary study of the onset of collective motion shows that one can go directly from a cohesive non-moving droplet to the incohesive phase and that in this region the onset of collective motion is roughly independent on $`\beta `$, yielding a vertical boundary.

More work should be devoted to the determination of the asymptotic phase diagrams mentioned above, as well as to a quantitative study of the onset of collective motion. Even in the simplest case of the non-cohesive Vicsek-type models, it can be second or first order, depending on the nature of the microscopic noise in the model . In both cases, we have started to uncover a rich interplay between collective motion, critical fluctuations, rotation modes, and shape dynamics in the transition region.

Not even mentioning the study of our model in three dimensions and its possible applications to particular real-world situations (e.g. biology, zoology, and robotics) the work presented here is probably only the beginning of the exploration of this new type of non-equilibrium systems.

## Appendix: Mean-field approach

To understand the interplay between alignment and body force, we simplify our system via a Hamiltonian model. The total energy is the sum of kinetic energy, two-body interaction energy, and a term related to the velocity alignment “force”. We thus write

$`={\displaystyle \frac{Nm}{2}}v_0^2`$ $`+`$ $`U(r_1,r_2,\mathrm{},r_N)`$
$``$ $`\alpha {\displaystyle \underset{i,ji}{\overset{N}{}}}\stackrel{}{v}_i.\stackrel{}{v}_j\stackrel{}{h}_0.{\displaystyle \underset{i}{\overset{N}{}}}\stackrel{}{v}_i.`$
where we have introduced an external field $`\stackrel{}{h}_0`$.

In the spirit of the mean-field approach, we assume that the fluctuations are negligible so that each component can be integrated independently and approximated by its averaged value. Thus the contibution from the two-body interaction becomes:

$`U(r_1,\mathrm{},r_N)`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}u\left(|\stackrel{}{r}_i\stackrel{}{r}_j|\right){\displaystyle \frac{N}{2}}U_0,`$
$`U_0`$ $`=`$ $`{\displaystyle _{R_\mathrm{c}}^+\mathrm{}}\beta \rho u(r)d^dr2a\rho \beta ,`$
where $`a`$ is a constant which depends on the potential. The alignment energy is computed with the coarse–grained velocity $`\stackrel{}{\phi }=\stackrel{}{v}`$ and the average number of neighbors $`\rho \pi r_0^2=\rho s`$ :

$$\alpha \underset{i,ji}{\overset{N}{}}\stackrel{}{v}_i.\stackrel{}{v}_j+\stackrel{}{h}_0.\underset{i}{\overset{N}{}}\stackrel{}{v}_i\left(\stackrel{}{h}_0+\alpha \rho s\stackrel{}{\phi }\right).\underset{i}{\overset{N}{}}\stackrel{}{v}_i.$$

In the following, we use $`\stackrel{}{h}_0+\alpha \rho s\stackrel{}{\phi }=\stackrel{}{h}`$.

Because of the separation of phase space variables, each term of the Hamiltonian contributes a factor to the partition function:
$``$ the two-body coupling :

$$\left(SbN\right)^N\mathrm{exp}\left(aN\rho \frac{\beta }{kT}\right),$$
(8)
where $`b=\frac{\pi }{2}R_c^2`$ is the hard–core surface and $`S`$ is the whole surface,
$``$ the alignment coupling :

$$\left[2\pi mv_0\mathrm{I}_0\left(\frac{hv_0}{kT}\right)\right]^N,$$
(9)
where $`\mathrm{I}_0`$ is a Bessel function. The free energy of the system is thus written:

$`{\displaystyle \frac{F}{NkT}}`$ $`=`$ $`\mathrm{ln}\left({\displaystyle \frac{mv_0}{\mathrm{}}}\right){\displaystyle \frac{mv_0^2}{2kT}}+1+\mathrm{ln}\left(\rho ^1b\right)`$
$`+`$ $`a\rho {\displaystyle \frac{\beta }{kT}}+\mathrm{ln}\mathrm{I}_0\left({\displaystyle \frac{hv_0}{kT}}\right)\rho s\phi ^2{\displaystyle \frac{\alpha }{2kT}}`$
Imposing to be at a minimum of the free energy, we find that the average velocity is the solution of a self–consistent equation:

$$\frac{\phi }{v_0}=\frac{\mathrm{I}_1}{\mathrm{I}_0}\left(\frac{v_0h}{kT}\right),$$
(10)
and that the critical point is defined by

$$sv_0^2\rho \alpha =2kT.$$
(11)
Collective motion thus emerges via a continous phase transition whose critical exponents are $`1/2`$ for the order parameter and $`1`$ for its susceptibility, independently of the cohesive force.

Furthermore, the necessary concavity of the free energy function implies a first order cohesion transition, since the second derivative of the free energy has two zeros. Their positions define the phase coexistence region. We determined the stability limit of the gas phase. As positions and velocities are decorrelated, we studied two cases : first without any motion, and then with collective motion. In the first case:

$$\beta _{\mathrm{GL}}=\frac{kT}{2a\rho \left(1b\rho \right)^2},$$
(12)
which means that the transition line does not depend on the alignment parameter (Fig. 16). Note also that there is no transition point in the zero–density limit of this model. We numerically solved the case with collective motion and found a stabilization of the “liquid” phase (Fig. 16). This is an effect of the assumption that velocity fluctuations are negligible. At the onset of motion, we expect that fluctuations of the macroscopic velocity diverge, leading to strong density fluctuations. This could explain the de–stabilization of the “liquid” phase (Fig. 15) and its stabilization in the mean field model 16, where such fluctuations are by definition absent.


text:
# Black-Hole Mass Measurements

## 1. Introduction

Knowledge of the mass of the central black hole in active galaxies is important for many studies of early structure formation and of the central engine and its evolution. The hope is that we will learn how the formation and growth of the black hole affects the formation, the evolution, and the characteristics of the galaxy by which it is hosted. The SDSS promises to provide valuable insight.

Black-hole mass determination methods fall in two distinct categories: primary and secondary methods. Their differences are important to recognize. With primary methods, the black-hole mass ($`M_{\mathrm{BH}}`$) is directly measured from gas and stars whose dynamics are dictated by the black hole. These methods are often very challenging or time consuming. In contrast, secondary methods are often more easily applied to large data sets, but only provide estimates of the mass by adopting approximations to the primary methods or by measuring parameters with which $`M_{\mathrm{BH}}`$ is known to correlate.

Due to space limitations, all the work done on this subject cannot be cited; my apologies in advance. Emphasis is placed on a few key papers providing calibration of the methods and assessments of the uncertainties.

## 2. Primary Mass Determination Methods

#### Stellar and Gas Kinematics.

To determine the virial black-hole mass of nearby quiescent galaxies requires high spatial resolution spectroscopy of the nuclear regions in order to measure the velocity dispersion of the stars (or gas) within the sphere of influence of the central black hole. This method can also be used for weakly active galactic nuclei (AGNs), such as low-luminosity AGNs, LINERs, and Seyfert 2s (Table 1). For more luminous AGNs (Type 1 sources) the strong nuclear emission washes out the stellar features in the spectrum. The scatter in the $`M_{\mathrm{BH}}`$$`\sigma _{}^{\mathrm{bulge}}`$ relationship, established for nearby quiescent galaxies (e.g., Tremaine et al. 2002), suggests that black-hole masses determined from stellar and gas kinematics are accurate to within a factor of about 2.

#### Megamasers.

Spectroscopy of water-vapor maser emission from circumnuclear disks in nearby AGNs reveals the kinematics of the disk and the location of the maser. Megamasers can yield highly accurate $`M_{\mathrm{BH}}`$ measurements, especially in highly inclined sources: for NGC 4258, which is inclined 83$`\pm `$4, the mass is determined to within a few per cent (Miyoshi et al. 1995). However, this method is only useful for selected ‘edge-on’ sources and so far relatively few objects are known to have megamasers (e.g., Greenhill et al. 2003).

#### Reverberation Mapping.

The reverberation-mapping technique is the best and most robust mass measurement method to apply to AGNs and quasars: utilizing the variability properties of the source the technique does not require high spatial resolution and is not affected by the nuclear glare. Peterson & Onken (2004, this volume) outline the basic principle of the reverberation mapping technique and discuss the importance of the zero-point calibration of reverberation masses and of the unknown, order unity, scale factor $`f`$.

For a given AGN the rms velocity width and time lag for different broad emission lines correlate such that higher ionization lines have larger widths and smaller lags (consistent with the same central mass for each object). This virial relationship for multiple lines is seen for all four AGNs for which this is testable (NGC 7469, NGC 3783, NGC 5548, 3C 390.3) and for the well-measured emission lines of Si iv $`\lambda `$1400, C iv $`\lambda `$1549, He ii $`\lambda `$1640, C iii\] $`\lambda `$1909, H $`\beta `$ $`\lambda `$4861, and He ii $`\lambda `$4686 (Peterson & Wandel 1999, 2000; Onken & Peterson 2002).

Notably, nearby AGNs with measurements of both reverberation mass and bulge velocity dispersion fall along the $`M_{\mathrm{BH}}`$ $`\sigma _{}`$ relationship established by quiescent galaxies and with a similar scatter (Ferrarese et al. 2001). This suggests that reverberation mapping masses are accurate to within a factor 2 to 3.

Reverberation mapping is fully applicable to distant quasars (Table 1), but is less practical as it is extremely time and resource consuming: luminous AGNs vary with smaller amplitudes and on longer time scales than lower-luminosity, nearby AGNs. The time scales are further increased by time dilation owing to the cosmological distances of quasars. Even for the 17 nearby quasars monitored so far, ten years of variability data are necessary to obtain reasonable results (Kaspi et al. 2000). Secondary methods are more practical for distant quasars.

## 3. Secondary Methods

Since the primary mass determination methods are either inapplicable or impractical for more luminous and more distant AGNs and quasars, several secondary methods have been adopted in the literature to estimate the central mass. These methods, summarized in Table 2, are either approximations to the reverberation mapping technique or rely on the empirical relationships between the black-hole mass, $`M_{\mathrm{BH}}`$, and the properties of the host galaxy bulge: velocity dispersion, $`\sigma _{}`$, or bulge luminosity, $`L_{\mathrm{bulge}}`$. The measurements of the bulge properties are most useful for AGNs at low redshift ($`z\stackrel{<}{_{}}1`$) where the host galaxy is easier to characterize. It is noteworthy that, even so, these methods can yield mass estimates that are currently quite uncertain, as explained below.

### 3.1. Scaling Relationships

Based on single-epoch spectroscopy the scaling relationships are approximations to the virial-mass measurements obtained from reverberation mapping. The method relies on the radius – luminosity relationship, established by reverberation mapping. For a photoionized BLR, its size $`R`$ should scale with the square root of the source luminosity, $`RL^{0.5}`$. Kaspi et al. (2000) found empirically that the size of the H $`\beta `$ emitting region scales with the continuum luminosity at 5100Å to the power of 0.7: $`R`$(H $`\beta `$) $`L_\lambda ^{0.7}`$(5100Å). While this difference is yet to be understood (work is in progress; e.g., Peterson & Onken, this volume), this means that $`M_{\mathrm{BH}}`$ can be estimated when we have measurements of the continuum luminosity and the emission-line width: $`M_{\mathrm{BH}}`$ $``$ FWHM2$`L_\lambda ^{0.7}`$. (Note, the slope adopted in the literature is not uniform but ranges between 0.5 and 0.7). The use of FWHM is a single-epoch approximation to the velocity dispersion of the line-emitting gas responding to continuum variations (Peterson & Onken 2004, this volume). Vestergaard (2004b, this volume) explains why these scaling relations are reasonable and why they can be applied to luminous, high-$`z`$ quasars.

Three different broad emission lines (H $`\beta `$, Mg ii, and C iv) have been adopted along with continuum luminosities (at 5100Å, 3000Å, and 1350Å, respectively) to estimate black-hole masses in large samples of quasars and AGNs \[H $`\beta `$, C iv: Vestergaard (2002; 2004b, this volume); Mg ii: McLure & Jarvis (2002); Jarvis & McLure (2004, this volume)\]. There are pros and cons to each method as outlined next. Using optical spectroscopy, the H $`\beta `$ method probes redshifts below 0.9, Mg ii probes redshifts between 0.3 and 2.3, while C iv is accessible at redshifts between 1.2 and 4.9. While Fe ii emission heavily contaminates the Mg ii line profile (e.g., Francis et al. 1991; Vestergaard & Wilkes 2001) and thus needs to be subtracted to allow measurements of FWHM(Mg ii), typically only strong Fe ii emission affects the H $`\beta `$ line width. C iv is only affected in the extreme wings by He ii$`\lambda `$1640, and occasionally by N iv\]$`\lambda `$1486, and (weak) Fe ii multiplets. However, the contribution to H $`\beta `$ from the narrow-line region needs to be subtracted. The only concern about C iv is the possible presence of outflowing high-ionization BLR gas, which may yield blue asymmetric broad profiles. This is seen very clearly in the case of narrow-line Seyfert 1 galaxies (e.g., Leighly 2000) for which the C iv method is not well suited. But very few quasars appear to display such asymmetric profiles, so emitting outflows are probably not a general concern for quasars (Vestergaard 2004). A notable assumption enters the Mg ii method, namely that Mg ii should be emitted co-spatially with H $`\beta `$ and should therefore have similar line widths (which seems to be the case; McLure & Jarvis 2002). Then FWHM(Mg ii) can be used as a direct surrogate for H $`\beta `$ in the scaling relation. Unfortunately, the only published data on Mg ii constrain the Mg ii lag (or the distance from which it is emitted) very poorly but the analysis by Clavel et al. (1991) and Dietrich & Kollatschny (1995) suggest that Mg ii is not necessarily emitted co-spatially with H $`\beta `$. The additional advantage of the available H $`\beta `$ and C iv scaling relations is that they are calibrated to the more accurate reverberation mapping masses of nearby AGNs (Wandel, Peterson, & Malkan 1999; Kaspi et al. 2000; Vestergaard 2002). The statistical scatter indicates 1$`\sigma `$ uncertainties of factors 2.5 to 3 (relative to the reverberation masses) for these scaling relations. Similar uncertainties are argued to be obtainable with Mg ii (McLure & Jarvis 2002). Nonetheless, individual mass estimates may be uncertain by as much as a factor 10. The power of scaling relations is in their application to large statistical AGN samples.

The uncertainties of these scaling relations are dominated by the scatter in the $`RL`$ relation, and improvements to this relation is thus desirable. Also, there is a need to understand the outliers in the scaling relations (Table 2). Dramatically improving the accuracy of the reverberation masses, on which the scaling relations rely, will require a dedicated space-based observatory, such as the proposed MIDEX mission Kronos (e.g., Peterson et al. 2004).

The radius $`R`$ can also be estimated with the reasonable assumption that the product of the ionization parameter and the electron density is roughly similar in all objects (the photoionization method). In this case the mass estimate is a function of $`L^{0.5}`$ FWHM2 and seems accurate to within a factor two of the reverberation masses (Wandel, Peterson, & Malkan 1999).

### 3.2. Via the $`𝑴_{\mathrm{𝐁𝐇}}\mathbf{}𝝈_{\mathbf{}}`$ Relationship

As the $`M_{\mathrm{BH}}\sigma _{}`$ relation for quiescent galaxies is relatively tight and AGNs also fall on this relationship, good measurements of the bulge velocity dispersion in AGN host galaxies can be used to estimate the black-hole mass (and with a similar uncertainty; Table 2). So far, only AGNs at $`z<0.06`$ have had $`\sigma _{}`$ measured because the Ca ii$`\lambda `$$`\lambda `$8498, 8542, 8662 absorption lines used for this measurement move into the atmospheric water vapor bands at $`z0.06`$ (Ferrarese et al. 2001).

Nelson (2000) suggested using FWHM(\[O iii\]) as a proxy for $`\sigma _{}`$ since the near-nuclear stellar velocity dispersion correlates with FWHM(\[O iii\]) for a sample of 75 Seyfert galaxies (Nelson & Whittle 1996). Boroson (2003) used data from the SDSS Early Data Release to assess and confirm that the 1$`\sigma `$ uncertainty of this method is a factor of 5. For this method to be more useful, both the scatter and the outliers need to be understood (Table 2). Some line asymmetries are suspected to be connected with outflows and can be prominent in radio sources where the narrow-line region and the radio source may interact.

Studies in the literature have used AGN host galaxy imaging and the Fundamental Plane for elliptical galaxies to measure the effective surface brightness, $`_e`$, and the effective radius, $`r_e`$, to derive first $`\sigma _{}`$ and then $`M_{\mathrm{BH}}`$. The reasoning is that especially nearby quasars and radio galaxies seem to reside in giant elliptical galaxies, which should fall on the Fundamental Plane. Unfortunately, this method potentially has large uncertainties because (a) $`_e`$ and $`r_e`$ are very hard to measure accurately in the presence of the bright nucleus of most AGNs, (b) the bulge/disk decomposition process is difficult even for nearby quiescent galaxies and is particularly challenging for AGNs, and (c) the method is subject to the uncertainties of both the Fundamental Plane and the $`M_{\mathrm{BH}}`$ $`\sigma _{}`$ relation. For radio galaxies, unaffected by nuclear glare, a factor 4 (or larger) uncertainty in estimating $`\sigma _{}`$ alone seems appropriate (Woo & Urry 2002). The combined uncertainty of the Fundamental Plane method is thus expected to exceed a factor 5.

### 3.3. Via the $`𝑴_{\mathrm{𝐁𝐇}}\mathbf{}𝑳_{\mathrm{𝐛𝐮𝐥𝐠𝐞}}`$ and Scaling Relationships

The black-hole mass of quiescent galaxies correlates, in addition to $`\sigma _{}`$, with the bulge luminosity, $`L_{\mathrm{bulge}}`$, although with larger scatter. For this reason, the $`M_{\mathrm{BH}}`$ $`L_{\mathrm{bulge}}`$ relation is generally not preferred for mass estimates of normal galaxies. This method has also proven difficult and prone to significant uncertainties owing to the nuclear glare affecting the bulge/disk decomposition process (e.g., Wandel 2002). Nonetheless, McLure & Dunlop (2001, 2002) argue that with careful $`R`$-band measurements of the AGN host galaxy, first $`L_{\mathrm{bulge}}`$ and then $`M_{\mathrm{BH}}`$ can be estimated to within a factor 3 to 4. Unfortunately, this method has not yet been calibrated to reverberation mapping masses, but rather to mass estimates based on scaling laws. This may not be a severe problem since the authors find their AGNs to have the same slope and scatter as the inactive galaxies with dynamical mass measurements to within the errors.

## 4. Future Efforts

Table 2 lists suggestions to how we may improve the mass estimation methods. In addition, it is important to understand their limitations as well as their efficacy. Since reverberation masses provide the anchor of the scaling relations a few comments thereon are in order. An understanding of the BLR structure and kinematics is required to determine the absolute zero-point of the reverberation masses (the ‘$`f`$’ factor). This is a non-trivial but an achievable task with a dedicated observatory like Kronos and existing advanced analysis techniques (Peterson et al. 2004; Peterson & Horne 2004). Also, direct comparisons of stellar dynamics and reverberation mapping need to be effected. Moreover, the odd objects need to be understood: for example, for 3C 390.3 the lower ionization lines appear to have shorter lags, in contrast to other AGNs. As this object is the only broad-line radio source which has been monitored, it needs to be established whether or not it belongs to a separate class of objects. Understanding this behavior is also important for our understanding of the typical BLR structure.

## References

Boroson, T. A. 2003, ApJ, 585, 647

Clavel, J., et al. 1991, in Variability of Active Galactic Nuclei, ed. H. R. Miller & P. J. Wiita (Cambridge: Cambridge University Press), 301

Dietrich, M. ,& Kollatschny, W. 1995, A&A, 303, 405

Dietrich, M., et al. 2002, ApJ, 581, 912

Ferrarese, L., et al. 2001, ApJ, 555, L79

Francis, P., et al. 1991, ApJ, 373, 465

Greenhill, L., et al. 2003, ApJ, 582, L11

Kaspi, S., et al. 2000, ApJ, 533, 631

Leighly, K. 2000, New Astronomy Review, 44, 395

McLure, R., & Dunlop, J. 2001, MNRAS, 327, 199

McLure, R., & Dunlop, J. 2002, MNRAS, 331, 795

McLure, R., & Jarvis, M. 2002, MNRAS, 337, 109

Nelson, C. 2000, ApJ, 544, L91

Nelson, C., & Whittle, M. 1996, ApJ, 465, 96

Onken, C., & Peterson, B. M. 2002, ApJ, 572, 746

Peterson, B. M., & Horne, K. 2004, Astron. Nach., in press

Peterson, B. M., Polidan, R. S., & Horne, K. 2004, Astron. Nach., in press

Peterson, B. M., & Wandel, A. 1999, ApJ, 521, L95

Peterson, B. M., & Wandel, A. 2000, ApJ, 540, L13

Tremaine, S., et al. 2002, ApJ, 574, 740

Vestergaard, M. 2002, ApJ, 571, 733

Vestergaard, M. 2004, ApJ, 600, in press

Vestergaard, M., & Wilkes, B. J., 2001, ApJS, 134, 1

Wandel, A. 2002, ApJ, 565, 762

Wandel, A., Peterson, B. M., & Malkan, M. A. 1999, ApJ, 526, 579

Woo, J.-H., & Urry, C. M. 2002, ApJ, 579, 530


text:
# Optical properties and Raman scattering of vanadium ladder compounds

## I Introduction

The vanadium ladder structures NaV2O5, CaV2O5, and MgV2O5 are low-dimensional compounds with unusual physical properties due to a strong coupling of charge, spin, and lattice degrees of freedom. NaV2O5 and CaV2O5 are iso-structural (space group $`D_{2h}^{13}`$ $`(Pmmn)`$) while MgV2O5 has the symmetry group $`Cmcm`$.The main structural element of all these compounds is a simple ladder containing V-O rungs and V-O legs providing a possibility for electrons to hop within a rung and between the rungs of one ladder, whereas the hopping matrix element between different ladders is considerably smaller. The highest occupied states are formed by $`d_{xy}`$ electrons of V mixed with oxygen $`2p`$ orbitals. CaV2O5 and MgV2O5 are both half-filled compounds with similar electronic properties, i.e. one $`d_{xy}`$ electron per V site. However, they show different magnetic behavior Korotin et al. (2000). The quarter-filled NaV2O5 where one $`d_{xy}`$ is shared by two V sites within a rung Smolinski et al. (1998), undergoes a spin-Peierls transition at $`T_c=35`$ K Isobe and Ueda (1996) accompanied by a charge ordering and a large lattice distortion. In this compound, phonons induce strong charge fluctuations near $`T_c`$ that lead to fluctuations in the spin-spin exchange $`J`$ Sherman et al. (1999). Therefore, investigation of lattice dynamics and electron (spin) - phonon coupling in ladder compounds can provide a clue to their properties.

Investigations of these compounds by optical spectroscopies - ellipsometry and Raman scattering - have shown many interesting features. The former gives the dielectric function, whereas the latter provides information about phonons and low-energy excitations and their coupling to light. For NaV2O5 the phonon Raman spectra demonstrate very large changes below the transition point Fischer et al. (1999) such as variations in the electronic background Sherman et al. (2001) and the appearance of intense new Raman active modes indicating a strong coupling of the lattice to the order parameter of the transition Fischer et al. (1999); Sherman et al. (1999). With the aim to understand similarities and differences in the optical properties of NaV2O5, CaV2O5, and MgV2O5 we perform ab initio investigations within density functional theory (DFT). We calculate the fully symmetric Ag phonon modes and their influence on the electronic structure and the dielectric function. Such phonon-induced modulations of the polarizability are a measure for electron-phonon coupling and the source of phonon Raman scattering of light. Although the vanadates are considered to be correlated materialsMazurenko et al. (2002), the results obtained in the framework of DFT show good agreement with experiment in terms of their lattice dynamics and optical properties. Since similar findings have been reported earlier for high temperature superconductorsAmbrosch-Draxl et al. (2002) such investigations are justified.

## II Computational details

All calculations are performed within DFT using the full-potential augmented plane waves + local orbitals (FP-APW+lo) E. Sjöstedt and Singh (2000) formalism implemented in the WIEN2k code Blaha et al. (2000). Exchange and correlation terms are described within the generalized gradient approximation (GGA) Perdew et al. (1996). All atomic positions of the three structures have been relaxed as a starting point for the investigation of the phonon modes within the frozen-phonon approach. We have displaced the atomic coordinates from their equilibrium positions, where for each degree of freedom two displacements in both, positive and negative direction, respectively, have been used. The corresponding forces established the dynamical matrix for the Ag phonons. The dielectric tensors are computed within the random phase approximation with the Kohn-Sham eigenstates taken as an approximation to the single particle states.

## III Dielectric function

Fig. 1 and Fig. 2 demonstrate the diagonal in-plane components of the dielectric function $`\epsilon _{xx}(\omega )`$ and $`\epsilon _{yy}(\omega )`$ of the three compounds investigated. The real part of $`\epsilon (\omega )`$ was obtained by the Kramers - Kronig relation from its imaginary part, where no broadening of the optical transitions was assumed. The strong absorption peak in the imaginary part of $`\epsilon _{xx}(\omega )`$ at $`\omega `$ close to 1 eV and another weak feature in Im $`\epsilon _{yy}(\omega )`$ at $`\omega `$ $``$ 1.25 eV are due to transitions between the bands originating from bonding and antibonding states within one rung separated by $``$ 2$`|t_{}|`$ in energy, where $`t_{}`$ is the in-rung hopping matrix element. Due to the dispersion of the bands the peak position is larger than 2$`|t_{}|`$. Our results correspond to $`t_{}=0.35`$ eV, in agreement with Ref. Smolinski et al., 1998. At the same time, a peak at 1.75 eV in $`\mathrm{Im}`$ $`\epsilon _{yy}(\omega )`$ arises from transitions within the legs. For NaV2O5 the calculated dielectric functions are in good agreement with available experimental data Presura et al. (2000). CaV2O5 demonstrates similar features, while for MgV2O5 the interpretation is more complicated.

## IV Ag phonon modes

The Ag phonon modes preserve the symmetry of the unit cell and in particular the inversion symmetry, and, therefore, are seen in the Raman spectra. Table 1 presents the calculated and measured eigenfrequencies of these modes for all three compounds. In all modes both, V and O displacements are considerably involved in the eigenvectors. The apical oxygen ions only contribute to the highest-frequency modes ($`\omega _{\mathrm{ph}}>900`$ cm-1), thereby vibrating along the $`z`$-axis. The theoretical eigenfrequencies of NaV2O5 and CaV2O5 for all but the lowest mode are in very good agreement with experiment, while the agreement for the MgV2O5 modes is worse. The eigenvectors of all three compounds correspond well to the assignment given in Ref. Popovic et al., 1999.

## V Phonon-induced changes in the dielectric function

Raman scattering arises due to phonon-induced changes in the band structure, and, correspondingly, in the dielectric function $`\delta \epsilon _{\alpha \alpha }(\omega )`$, where $`\alpha `$ are the Cartesian indices. The intensity of the one-phonon scattering $`I_R`$ for the incident and scattered light parallel to the $`\alpha `$ axis can be written as

$$I_R\left|1|\frac{\epsilon _{\alpha \alpha }}{Q}\widehat{Q}|0\right|^2.$$
(1)
Here $`Q`$ is the coordinate of the phonon mode, and the matrix element represents transitions between the states $`|0`$ and $`1|`$, which are the phononless ground state and the one-phonon state, respectively.

Fig. 3 represents the ion displacements in the V–O2 stretching mode which is taken as an example to demonstrate the influence of the phonons on the dielectric function. Figs. 4 and 5 show the changes in the $`xx`$ and $`yy`$ components of the dielectric function of NaV2O5 and CaV2O5, respectively. The black lines correspond to the situation of the undisplaced positions, whereas the grey lines reflect the distorted lattice with the displacement of the O2 ion along the $`x`$-axis taken to be 0.02 Å, which is of the order of magnitude of the zero-point vibrations. The V ion is shifted correspondingly according to the eigenvector of the mode. As it is seen in the figure, this mode causes big changes in $`\epsilon _{xx}(\omega )`$ at $`\omega `$ close to the resonance at 1 eV which corresponds to the intra-rung transitions due to a modulation of the $`t_{}`$ hopping within a rung. A general feature of the phonon-induced changes in the dielectric functions is that all the peaks are considerably shifted. This fact implies large values for the first derivatives of the dielectric functions, which govern the Raman scattering intensities. We therefore predict strong resonance effects in phonon Raman scattering close to these peak positions. E.g. the 1 eV region of NaV2O5 should exhibit resonances in the $`xx`$ component. For CaV2O5 we expect even more resonances since the peak structure of the dielectric function is more pronounced. Due to our findings we suggest Raman scattering experiments to be performed below 2 eV, the energy region where we expect very high Raman intensities. At the same time, we can explain the experimentally observed intensity enhancement for NaV2O5 when changing the frequency of the laser beam from 1.95 to 2.4 eV in $`xx`$ polarizationFischer et al. (1999), since the latter frequency is closer to the peak in Re $`\epsilon _{xx}(\omega )`$ at $``$3 eV as can be seen in the upper panel of Fig.1.

## VI Conclusions

We have performed ab initio calculations for the optical properties and Ag phonon modes in the ladder compounds NaV2O5, CaV2O5, and MgV2O5. For NaV2O5 the calculated $`\epsilon (\omega )`$ compares well with available experimental data. For NaV2O5 and CaV2O5 the phonon frequencies are in a very good agreement with the experiment, while for MgV2O5 the agreement is not as good. Due to the strong influence of the ion displacement on the dielectric function we predict pronounced resonance effects in phonon Raman scattering in the energy range below 2 eV.

Acknowledgment

We acknowledge support from the Austrian Science Fund (project P15520), and the EU (RTN network EXCITING, contract HPRN-CT-2002-00317).