Statistical mechanics and numerical modelling of geophysical fluid dynamics

(1)

UvA-DARE is a service provided by the library of the University of Amsterdam (http

s

://dare.uva.nl)

Dubinkina, S.B.

Publication date

2010

Link to publication

Citation for published version (APA):

Dubinkina, S. B. (2010). Statistical mechanics and numerical modelling of geophysical fluid

dynamics.

General rights

It is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s)

and/or copyright holder(s), other than for strictly personal, individual use, unless the work is under an open

content license (like Creative Commons).

Disclaimer/Complaints regulations

If you believe that digital publication of certain material infringes any of your rights or (privacy) interests, please

let the Library know, stating your reasons. In case of a legitimate complaint, the Library will make the material

inaccessible and/or remove it from the website. Please Ask the Library: https://uba.uva.nl/en/contact, or a letter

to: Library of the University of Amsterdam, Secretariat, Singel 425, 1012 WP Amsterdam, The Netherlands. You

will be contacted as soon as possible.

(2)

Chapter 4 A thermostat closure for point

vortices

4.1 Background

Invis id uid models are natural in a number of appli ation areas, su h as

atmosphereand o eans ien e, wheretheReynoldsnumbers areso largeas to

be ee tively innite. These ows are hara terized by onservation of total

energy,the as adeofvorti itytoeverners ales,andsensitivedependen eon

initial onditions [74 ℄.

For the numeri al simulation of su h ows, the la k of a vis ous diusion

lengths alepresentsthe hallengethat duetothevorti ity as ade,anydire t

dis retizationoftheequationsofmotionmusteventuallybe omeunderresolved,

as vorti ity istransported to s alesbelow thegridresolution. Ittherefore

be- omes ne essary to lose thenumeri almodel by some means. Any nite

nu-meri aldis retizationimpliesa losureofsomekind,whetherexpli itlymodeled

or impliedbythedis retization [21 ℄.

The most ommon approa h is the introdu tion of arti ial vis osity,

ei-therthroughmodi ationoftheuidequationstoin lude(hyper-)vis osity,or

throughtheuseofstabilizeddis retizations,forwhi hthevis oustermsappear

in a modied equation analysis [38 ℄. In either ase the vis ous length s ale

must beon theorder of thegrid resolution to be ee tive. Onedisadvantage

withavis ous losuremodelisthatitpre ludesanups ale as adeofvorti ity,

thereby suppressinghydrodynami instabilities in geophysi al ows,espe ially

in three dimensions. Alternatively, methods an be onstru tedthat preserve

the dis rete totalenergy exa tly. However,this isa hieved via a nonphysi al

re-inje tion oftheenergyfrom sub-grids alevorti ityat thelarges ales[63℄.

Aproper losuremodelshoulddistinguishbetweenresolvedandunresolved

s alesanda ountfortheex hangebetweenthese. Inthis hapterwe onsidera

simpletwo-s alepointvortexow, onsistingofasmallnumberofvorti eswith

large ir ulation and a very largenumber of vorti es with a mu h smaller

ir- ulation. Weseekasimplied omputationalmodelfortheaggregatebehavior

(3)

ble, in whi h a system of parti les is in thermal equilibrium with a reservoir.

This point of view hasbeenexploited by Bühler [9℄in a numeri al/statisti al

investigation of the work of Onsager [65 ℄, and our goal here is to reprodu e

theresultsof[9℄withoutexpli itlya ountingfortheindividualmotionsofthe

reservoirofsmallvorti es.

A thermostat is a tool used in mole ular dynami s to model a system in

thermalequilibriumwithareservoir;su hthermostatsmaybeeithersto hasti

(e.g. Langevin dynami s) or deterministi . In a Langevin dynami s

simula-tion a sto hasti perturbation is introdu ed in the for e eld together with a

dissipative term; these termsare maintained in balan e so as to preserve the

anoni al ensemble. With a dynami al(or deterministi )thermostat, by

on-trast, the system is augmented by a few degrees of freedom that model the

ex hange with thereservoir. The goalof thermostating isto for e thesystem

tosamplethe anoni alequilibriumdistributionatagiventemperatureby

on-tinuallyperturbingit. A benetof thedynami almodels isthat itis possible

to onservestru ture(e.g. Hamiltonianstru ture)in theaugmenteddynami s.

A motivation for using this approa h is that if the perturbation is small, the

dynami swillstill orrespondtophysi aldynami s(in onta twithareservoir)

onanintermediatetimes ale. Someexamplesofimportantdeterministi

ther-mostatingmethodsaretheNosémethod[61 ,62 ℄,whi hpreservesHamiltonian

stru ture at the expense of a ontinuous res aling of time, the Nosé-Hoover

method [62 ,37℄whi hre oversthelineartimebutloses anoni alHamiltonian

stru ture, the Nosé-Poin aré method of Bond, Laird & Leimkuhler [6 ℄ whi h

is anoni allyHamiltonian, and a generalization ofthe Nosé-Hoover approa h

forHamiltoniansystemswithPoissionstru ture[8℄. Deterministi thermostats

havealsobeen oupledwithLangevinmodelsin[44 ℄forexample.

Inthis hapterweproposetheuseofathermostat tomodeltheunresolved

vorti ityanditsex hangein asimplepointvortexmodel. Wewillmodel both

aninnitereservoirasin lassi althermodynami s,andanitereservoirashas

been used in theexperiments of [9 ℄. Statisti s of the thermostated dynami s

will be ompared with the results of Bühler [9 ℄. In 4.2 we make use of a

generalized thermostat whi h an be used to for e a Hamiltonian system to

sampleageneral lassofequilibriumdistributions. Thepointvortexmodeland

its statisti al me hani s is reviewed in 4.3. In 4.4 we present thedetails of

thethermostatednumeri almethods onsidered,in ludingthemodelsfornite

and innitereservoirs. Finally, in 4.5the numeri al s hemes are veried by

omparisonwithresultsfromtheliterature.

4.2 Generalized thermostats

Consider anopensubset

D ⊂ R

d

andadeterministi dierentialequation

˙

(4)

4.2. Generalized thermostats

79

A probability distribution

ρ(X, t)

∈ D × R → R

,

ρ

≥ 0

, on

D

is transported

under theve toreld

f

a ordingtothe ontinuityequation

∂

∂t

ρ(X, t) +

∇ · ρ(X, t)f(X) = 0.

(4.2)

This ontinuity equation implies that

R

D

ρ dX = 1

for all

t > 0

if this holds at

t = 0

. Anequilibrium distribution is a stationary solutionof(4.2). In this

hapterwewillbe on ernedprimarilywith systemsoftheform

˙

X

= J(X)

∇H(X),

X(t)

∈ D, J

T

₌

−J, H : D → R.

(4.3)

Thefun tion

H

isa rstintegral of(4.3),typi allytheenergy. If

J

is

indepen-dent of

X

, then this denes a (generalized) Hamiltonian system. Otherwise,

one mustalso show that

J(X)

satises theJa obi identity, in whi h asethe

systemisPoisson. Wemaketheweakerassumptionthattheve toreldonthe

rightside of (4.3) is divergen e-free, i.e.

∇ · f(X) ≡ 0

, so that the transport

equation(4.2)simpliestotheLiouvilleequation

d

dt

ρ(X, t) =

∂

∂t

ρ(X, t) + f (X)

· ∇ρ(X, t) = 0.

(4.4)

Anequilibriumdistribution

∂ρ

∂t

≡ 0

mustsatisfy

f (X)

_{· ∇ρ(X) ≡ 0.}

Note that any fun tion

ρ(X) = ρ(H(X))

that depends on

X

through a rst

integral is an equilibrium distribution. If (4.1) has additional rst integrals

I

2 (X)

,...,

I

p

(X)

,thenanydistribution

ρ(H, I

2 , . . . , I

p

)

isalsoanequilibrium

distribution. The ensemble average of a fun tion

F (X)

with respe t to the

equilibrium distribution

ρ(X)

is

hF i ≡

Z

D

F (X)ρ(X) dX.

Giventheiramplesupply,thedegreetowhi hagivenequilibriumdistributionis

meaningfullargelydependsonwhetherthesolutiontothedierentialequation

is ergodi in that distribution su h that thelong timeaverageof anyfun tion

F (X(t))

ofthesolution

F

_{≡ lim}

T →∞

1 T

Z

T

0 F (X(t)) dt,

onvergesto theensembleaveragein thedistribution,i.e. satises

F =

_{hF i}

(5)

Themi ro anoni alensemble[39 ℄appliestoanisolatedsystemat onstant

energy,andisthesingularmeasureontheenergylevelset ontainingtheinitial

ondition

ρ

ν

∝ δ(H(X) − H

0 ),

(4.5)

where

H

0 = H(X(0))

. Thisensembleisappropriateforanumeri alsimulation

withanenergy onservingdis retization.

A system in onta t with a large reservoir does not onserve energy, but

rather ex hanges it with the reservoir. If it is in thermal equilibrium with a

reservoir of statisti al temperature

β

−1

, then theappropriate ensemble is the

anoni alensemble[39 ℄ withGibbsmeasure

ρ(X) = N

−1

_{exp [}

−βH(X)] ,

(4.6)

where

N =

R

D

exp [

−βH(X)] dX

. It is lear, however, that a single

solu-tion of the system(4.3) will not beergodi in the Gibbs measure, sin e with

probability one it will sample the onstant energy surfa e ontaining the

ini-tial ondition,whereas(4.6)assignsnonzeroprobabilityto allenergysurfa es.

Instead, to modela systemin thermal equilibrium with a reservoir, one must

devise a method whose dynami s samples phase spa e with probability given

by the anoni aldistribution (4.6). Thedevelopment ofmethods that dojust

this onstitutesana tiveeld of resear h. Anumber of te hniques have been

developedforsamplinginagivendistribution,in ludingMonteCarlos hemes,

whi h generate random ongurations or traje toriesa ording to the hosen

distribution; Langevin thermostats, in whi h the original system of ordinary

dierential equations is augmented by sto hasti for ing and generalized

dis-sipation terms; and deterministi thermostats, in whi h the reservoir itself is

modelled using a small number of additional degrees of freedom. The latter

approa heshavetheadvantage thattheygenerateplausible solutionbehavior,

and an beusedto ompute orrelations.

In the next two se tions we des ribe generalized Langevin dynami s and

generalizedsto hasti Bulga -Kusnezovthermostatsforsamplinginawide lass

ofequilibrium distributionsforHamiltoniansystems.

4.2.1 Langevin thermostat

If one integrates(4.3)numeri allyusing a symple ti integrator, the

Hamilto-nian willtypi allybewell- onserved. Asa result,the solutionwill notsample

phasespa ewiththemeasure(4.6)above,but insteadwillstayneartheinitial

energy level set (at best sampling

ρ

ν

). For some appli ations it is desirable

to onstru ta perturbeddynami alsystemthatdoessample

ρ

whileretaining

something ofthedynami albehaviorof (4.3). Inthis wayone an onstru ta

plausible (representative) behavior of thesystem ifit were ex hanging energy

withthereservoir a ordingto

ρ

.

One approa h to sample a given equilibrium distribution augments (4.3)

with arefullytunednoise anddissipationterms:

˙

(6)

4.2. Generalized thermostats

81

where

g(X) :

D → R

d

,

Λ(X)

∈ R

d×d

isamatrixvaluedfun tion,and

w(t)

isa

ve torWiener pro ess,i.e. the

w

i

(t)

,

i = 1, . . . , d

, ares alar Gaussianrandom

variableswithmeanzeroandin rements

w

i

(t)

−w

i

(s)

∼ N (0, t−s)

. Phasespa e

densities are transported by the ow of (4.7) a ording to the Fokker-Plan k

equation(see,forexample[66 ℄)

∂

∂t

ρ(X, t) =

−∇ · ρ(X, t)(f(X) + g(X)) +

1

2 ∇ · ∇ · ρ(X, t)Λ(X)Λ

T

_{(X), (4.8)}

where

g(X)

mustbedeterminedsu hthatthedesiredequilibriumdistributionis

a stationarysolutionof(4.8). If

ρ

dependson

X

onlythroughitsHamiltonian

ρ(X) = ρ(H(X))

, then the Hamiltonian dynami s drops out of the

Fokker-Plan kequation,sin e

∇ · ρ(H(X))J∇H(X) = ρ∇ · J∇H − ρ∇H · J∇H = 0

bythedivergen efreenatureoftheHamiltonianowand onservationofenergy.

TheFokker-Plan kequationredu esto

∂

∂t

ρ(X, t) =

−∇ · ρ(X, t)g(X) +

1

2 ∇ · ∇ · ρ(X, t)Λ(X)Λ

T

_(X).

Therighthand sideiszeroif

−ρ(X)g(X) +

1 ₂

∇ · ρ(X)Λ(X)Λ

T

(X) = 0.

For the anoni al equilibrium distribution (4.6), we an solve this for

g(X)

(suppressingthedependen eon

X

inthenotation):

ρg(X) =

1

2 ρ

∇ · ΛΛ

T

−

1 ₂

βρΛΛ

T

∇H = ρ

µ 1

₂

∇ · ΛΛ

T

−

1 ₂

βΛΛ

T

∇H

¶

.

For

Λ(X)

onstant,thisbe omes

g(X) =

₋

1

2 βΛΛ

T

∇H(X).

Hen e,theLangevindynami sis

˙

X

= J

_{∇H(X) −}

β

2 ΛΛ

T

∇H(X) + Λ ˙w.

(4.9)

If the ow map is in addition ergodi with respe t to

ρ

, then the generalized

Langevin dynami s (4.7) an be used to sample the anoni al distribution at

inversetemperature

β

.

Remark.

For

Λ = Λ(X)

lo allydenedthenoiseismultipli ative,onemust

(7)

4.2.2 A generalized Bulgac-Kusnezov method

Thefollowingapproa h generalizestheBulga -Kusnezovmethod[8℄andoers

additional exibility. The method has been proposed for anoni al sampling

in the mole ular dynami s setting in [43℄; here wetreat an arbitrary smooth

ensembleand applyitto theuidvortexmodel. Weintrodu ea newvariable

ζ

∈ R

and fun tions

s(X, ζ) :

D × R → R

d

and

r(X, ζ) :

D × R → R

andform

the oupledsystem

˙

X

= J

∇H(X) + s(X, ζ),

(4.10)

˙ζ = r(X, ζ).

(4.11)

We ask that the following extended measure beinvariant under the Liouville

equation:

˜

ρ(X, ζ)

_{∝ exp [−βF (X) − αG(ζ)]}

(4.12)

for

F

and

G

appropriatelydened fun tions. Inthe aseof (4.6)wewilltake

F

_{≡ H}

,butwe onsiderthismoregeneralformulationfornow. Notethatafter

integrationover

ζ

,thismeasureisoftheform(4.6). Thestationarity ondition

forthetransportequation(4.2)is

∇ · ˜

ρ (f + s) + ∂

ζ

(˜

ρr) = 0,

with

f = J

∇H(X).

Some al ulations give

0 = (f + s)

_{· ∇˜}

ρ + ˜

ρ

_{∇ · (f + s) + r}

∂

∂ζ

ρ + ˜

˜

ρ

∂

∂ζ

r

=

_{−β ˜}

ρ

_{∇F · (f + s) + ˜}

ρ

_{∇ · (f + s) − α˜}

ρr

∂

∂ζ

G + ˜

ρ

∂

∂ζ

r

= ˜

ρ

µ

−β∇F · (f + s) + ∇ · s − αr

_∂ζ

∂

G +

∂

∂ζ

r

¶

,

wherethedivergen e-freedomoftheHamiltonianve toreldisusedinthelast

inequality.

Nextwemakesome simplifyingassumptions. Firstweassumethe

thermo-statvariable

ζ

tobenormallydistributed,taking

G(ζ) = ζ

2 _/2

. Wealsoassume

that

r

depends only on

X

, i.e.

r(X, ζ) = r(X)

. The stationarity ondition

onsequentlyredu esto

0 =

−β∇F · (f + s) + ∇ · s − αrζ.

Wewishtousethisrelationtodene

r

. Notethat

ζr(X) =

1 α

(

∇ · s − β∇F · (f + s)) .

(4.13)

Sin e

ζ

maybezero,ea htermontherightshouldeithervanishorhavepre isely

(8)

4.3. Statistical mechanics of point vortices

83

havefun tionaldependen e viatheHamiltonian. Ifweassumethat

F (H(X))

,

thentheskew-symmetry of

J

implies

∇F · f = F

′

(H(X))

_{∇H · J∇H ≡ 0.}

Ifadditionallyweassume

s(X, ζ)

tobelinearin

ζ

,i.e.

s(X, ζ) = s

1 (X)ζ,

s

1 (X)

∈ R

d

,

thenwendthat

r(X) =

1 α

(

∇ · s

1 (X)

− β∇F · s

1 (X))

(4.14)

isa solutionof(4.13).

Spe i hoi es of the fun tions

F (X)

and

s

1 (X)

will be treated in

Se -tion4.4.

Ingeneral,thethermostateddynami ssodened willnotbeergodi inthe

invariantmeasure(4.12). Toimproveergodi ity,aLangevintermmaybeadded

to (4.11). Seealso[44 ℄.

˙

X

= J

_{∇H(X) + s}

1 (X)ζ,

(4.15)

˙ζ = r(X) −

αλ

₂

2 ζ + λ ˙

w.

(4.16)

Sin ethenoiseentersthrough

ζ

,itinuen es

X(t)

onlyafterintegration,soits

ee t onthedynami sissmoothed.

Remark.

In the important spe ial ase

F (X) = F (H(X))

, ifwe hoose

s

1

su h that

∇ · s

1 ≡ 0

, thenthe system(4.15)-(4.16) an be astin theform

ofageneralizedLangevinthermostat(4.9)asdis ussedinthepreviousse tion.

Dene theaugmentedsystem

˜

X

=

Ã

X

ζ

!

,

J( ˜

˜

X) =

"

1 F

H

J

β

α

s

1 (X)

−

β

α

s

1 (X)

T

₀

#

,

˜

H(X, ζ) = F (H(X)) +

α

2β

ζ

2 _.

(4.17)

Then(4.15)-(4.16) with(4.14)takestheform

d

dt

X

˜

= ˜

J

∇ ˜

H

−

α

2 ΛΛ

T

∇ ˜

H( ˜

X) + Λ ˙

w

(4.18)

with

Λ =

¡

0 λ

¢

.

4.3 Statistical mechanics of point vortices

Themotionof

M

pointvorti eswith ir ulationstrengths

Γ

i

∈ R

,

i = 1, . . . , M

,

andpositions

x

i

(t)

∈ R

2

isgivenbytheHamiltoniansystem

Γ

i

˙x

i

= K

∂H

∂x

i

(9)

where

K = (

0 1

−1 0

)

,andtheHamiltonian

H =

₋

1 4π

X

i<j

Γ

i

Γ

j

ln(

|x

i

− x

j

|

2 )

representsthekineti energy.

Ifthereare

Γ

i

withbothpositiveandnegative ir ulations,thenthemotion

of point vorti es is unbounded on theplane. A bounded ow an be ensured

by imposing periodi ity, whi h altersthe Green'sfun tionin theHamiltonian

[86 ℄. Alternatively,owonadis ofradius

R

anmemodeledbydeningaset

ofimage vorti es

Γ

′

i

=

−Γ

i

,

x

′

i

= x

i

R

2 |x

i

|

2 ,

i = 1, . . . , M

whi hensurethatthevelo ityeld observedbyanypointvortexistangentto

the wall. Inthe dis model, whi h weadopt in this hapter,the Hamiltonian

hasthreetermsdueto: theoriginalpairpotential,theself-intera tion,andthe

intera tion termsofea h vortexwiththeimagesoftheothers:

H =

−

_4π

1 X

i<j

Γ

i

Γ

j

ln(

|x

i

− x

j

|

2 ) +

1 4π

X

i

Γ

2 i

ln(R

2 − |x

i

|

2 )+

1 4π

X

i<j

Γ

i

Γ

j

ln(R

4 − 2R

2 x

i

· x

j

+

|x

i

|

2 |x

j

|

2 ).

(4.20)

To astthesystem(4.19)in theform (4.3),wedene

X

= (x

T

1 ,

· · · , x

T

M

)

T

,

H = H(X)

, and

J =







Γ

−1

1 K

. . .

Γ

−1

M

K







.

Besides thekineti energy, thepointvortexow onthe dis onservesthe

totalangularmomentum,denedas

ˆ

M =

1 2π

X

i

Γ

i

|x

i

|

2 .

(4.21)

Ingeneraltherewillbeanex hangeofmomentumbetweenthestrongvorti es

andthereservoir. Howeveronaveragewewouldexpe ttheangularmomentum

of bothstrongand weak vortex setsto beapproximately onstant. Infa t, it

wouldbestraightforwardtomodeltheex hangeofangularmomentumusingthe

thermostataswell. Thiswouldrequireknowledgeofthevarian eoftheangular

momentumofthereservoir.Inthis hapterweassumethemomentumex hange

(10)

4.3. Statistical mechanics of point vortices

85

signi ant drift in angular momentum. To orre t this, one ould onstru t

a proje tion of the noise term

Λ(X)

onto the angular momentum manifold.

However,this omesatthe ostofmultipli ativenoise.

Thephase spa e of the point vortex ow onsists of thedire t produ t of

M

opiesofthedomain. Ifthedomainisbounded,so isthephasespa e. The

energy

H

isunboundedonthephasespa ehowever: as

x

i

→ x

j

,thelogarithm

tendsto

−∞

;if

Γ

i

and

Γ

j

arelike-signed,

H

→ +∞

,ifoppositely-signed,

H

→

−∞

. Inparti ular, ifa parti le ollides withthewall,

H

→ −∞

. Asnotedby

Onsager[65 ℄,ifwedene

Ω(E)

tobethemeasureofthesetof ongurationsin

phasespa eforwhi h

H

∈ [E, E +dE]

,thenwemusthave

lim

E→±∞

Ω(E) = 0

.

In other words, sin e the phase spa e is bounded, the measure of available

phasespa emust eventuallyde reaseas anin reasingfun tionofenergy. The

situationisin ontrasttoother

n

-bodyproblemsen ounteredin hemistryand

astronomy, where the positive denitekineti energyterms an a ommodate

anyamountofenergy, andthemeasure ofavailablephasespa eisa monotone

in reasingfun tionofenergy.

Consequently, the mi ro anoni al entropy

S(E) = ln Ω(E)

must attain a

maximum for some

E

∗

. The mi ro anoni al temperature is dened to be

T

−1

ν

=

dE

d

ln S(E)

. Astemperaturevaries from

−∞

to

∞

alongthe realline,

theasso iated energystatesvary from

E > E

∗

to

E < E

∗

through innity. In

other wordstemperature passes through zerovia a ollision. This has

impor-tant onsequen esforthermostating,sin eitimpliesthatthatthevorti es will

ollapsetothewallifthetemperature hangessign ontinuously.

Re all that extreme values of the energy

H

are asso iated with lose

ap-proa hes between vorti es or image-vorti es. As noted by Bühler [9 ℄, for a

homogeneous system with

Γ

i

= Γ

, the energy largely governs the dynami s,

sin e ollisionshavetoo urroughlyat onstantenergy. Thesituationismore

interestinginaheterogeneoussystemwithvorti esofgreatlydieringstrength.

Onsager predi ted that for su h systems, extreme values of energy would

in- rease the probability of lustering of like-signed or opposite-signed vorti es,

with a preferen e for the strongest ones, su h that mostof the energy would

reside in a fewdegrees offreedom. As a result,thesmallvorti es would roam

aimlessly about, not developing into oherent stru tures, but exhibiting large

entropy.

Bühler dis ussesOnsager's ideas in the ontext of the anoni al ensemble

appliedtothestrongvorti es,whi h onstituteasystemin`thermal'equilibrium

with the reservoir of weak vorti es. He veries Onsager's predi tions using

numeri alexperimentswithasystemof100pointvorti es,fourhavingstrength

±10π

andthe rest having strength

±2π

. In ea h group, half thevorti es had

positive ir ulationandhalfnegative. Experimentswere arriedoutforextreme

positive, neutral and extreme negative inverse statisti al temperatures

β

ν

=

T

−1

ν

inthemi ro anoni alsense. Inea h asethestrongvorti eshadthesame

(nearlysteadystate)initial onguration,sothedieren esinenergywereonly

(11)

and opposite signedstrongvorti es,distan e fromthe wall,and energyin the

strongvorti es.Bühlerdistinguishesbetweenatheoreti alinnitereservoirand

the nitereservoir omprisedofthe 96weakvorti es. In theinnitereservoir

ase,the anoni alprobabilitymeasureonlyexistsforaniteintervalofinverse

temperature

β

,whoseboundary orrespondsto ollisions. Thissituationisdue

to the availability of an innite amount of energy in the reservoir, and has

impli ations for thermostating in the anoni al ensemble. Spe i ally, if one

thermostats in the anoni al ensemble and in reases

β

beyond its admissible

range, the vorti es will ollapse onto the boundary. Bühler also points out

that the onta t with a nite reservoir will suppress this ollapse, allowing

thermostating at all temperatures. This is be ause there is a nite amount

of energy in the nite reservoir, and this ee tively bounds the loseness of

approa h of anytwo vorti es from below. Theprobability ofa loseapproa h

be omesverysmall. Theprobabilityof

H = E

fora systemin onta twitha

nite reservoirde ayslike

exp(

−γE

2 ₎

with

E

forsome

γ > 0

.

4.4 A thermostated integrator for point vortices

Our goal is to apply the generalized Bulga -Kusnezov thermostat from

Se -tion 4.2.2to thepointvortexow ofSe tion 4.3. Inthis se tion wellin the

details ofthemethod. First,in Se tion4.4.1wespe ifytwoequilibrium

distri-butions orrespondingtothe aseswherethereservoirofsmalls alevorti ityis

nite orinnite. InSe tion4.4.2wedene athermostat fun tion

s

1

su hthat

thegeneralizedBulga -KusnezovthermostatisaLangevinthermostat. We

de-s ribethenumeri almethodusedto integratethemodeladaptivelyintimein

Se tion 4.4.3andthemeansof omputingthetemperatureinSe tion4.4.4.

4.4.1 Infinite and finite reservoir ensembles

As dis ussed in [9 ℄ the behavior of a thermostated point vortex system an

vary onsiderably depending on whether the reservoir is nite or innite. In

the aseof an innitereservoir, as thetemperatureof the reservoir is pushed

toward zero, the subsystem may draw an arbitrarily large amount of energy

from thereservoir,leading to ollisionsbetweenindividualvorti es orwiththe

wall. For a nitereservoir,there is a limitedamountof energyavailable su h

thata ollisionmayonlyo ur ifa ollisionwithoppositeenergyo ursatthe

same time,and thisisimprobable. Spe i ally, inthe aseofa nitereservoir

with normally distributed reservoir energy, the equilibrium distribution takes

theform

˜

(12)

4.4. A thermostated integrator for point vortices

87

Forthegeneralizedthermostat (4.10)(4.11)we anmodelbothniteand

in-nitereservoirs. For anitereservoir wetake

F (X) = H(X) +

γ

β

H(X)

2 _,

r(X) =

1 α

(

∇ · s

1 (X)

− (β + 2γH(X))∇H · s

1 (X)) ,

(4.22)

andforaninnitereservoir

γ

≡ 0

in theexpressionsabove.

4.4.2 Choice of s

1

Wemakethefollowing hoi eforthefun tion

s

1

in (4.14):

s

1 (X) =

−







K x

1 |

x

1 |

. . .

K x

M

|

x

M

|







.

(4.23)

Theowoftheve tor eld

s

1

preservesthedistan eofea hpointvortexfrom

the enter ofthedomain. Consequentlythethermostated system(4.15)-(4.16)

preservestheangularmomentum(4.21).

Furthermore, this hoi eof

s

1

isdivergen e-free:

∇ · s

1 (X)

≡ 0,

implying that the thermostated dynami s is a generalized Langevin system

(4.18), and that the integral

H

˜

in (4.17) with

F (H(X))

from (4.22) is

pre-servedin thelimit

λ

→ 0

ofzeronoiseanddissipation.

4.4.3 Implementation details

Inournumeri alimplementation,timesteppingwas doneusing asplitting

ap-proa h. We solved alternately the deterministi thermostat system and the

sto hasti equation for the thermostat variable. The deterministi system is

solvedwiththeimpli itmidpointrule

X

n+1

_{− X}

n

τ

= J

∇H( ˆ

X

)

− s

1 ( ˆ

X

)ˆ

ζ,

(4.24)

ζ

n+1

_{− ζ}

n

τ

= r( ˆ

X),

(4.25)

where

X

ˆ

= (X

n+1

_{+ X}

n

_)/2

and

ζ = (ζ

ˆ

n+1

_{+ ζ}

n

_)/2

.

Theremainingve toreldisanOrnstein-Uhlenbe kequation

(13)

withexa t solution

ζ

n+1

= exp[

_−ετ]

Ã

z

n

+ λ

r

e

2ετ

_{− 1}

2ε

∆w

!

,

(4.26)

where

ε = αλ

2 _/2

and

∆w

∼ N (0, 1)

. A fulltimestepofsize

∆t

is onstru ted

by solving (4.26) with

τ = ∆t/2

omposed with (4.24)(4.25) with

τ = ∆t

omposedwithase ond solutionof(4.26),

τ = ∆t/2

.

Duringa loseapproa hoftwovorti es,equivalentlywhenthestrongvortex

energyislargeinmagnitude,a ura yandstability onsiderationsmotivatethe

useofanadaptivetime-steppingstrategy. Givenastepsize

∆t

n

inthe

n

thtime

step,thesubsequenttimestepisfoundbysolving

∆t

n

∆t

n+1

= e(X

n

)

2 ∆s

2 .

(4.27)

Here,

∆s

isauniformtimestepunderthetimetransformation

t = e

·s

,and

e

isa

monitorfun tionthatmeasuresthestinessofthelo alsolution. This

adaptiv-ityapproa h is expli it andtime-reversible whenever thenumeri alintegrator

issymmetri . For ourexperimentsweuse

e(x) = min

i6=j

|x

i

− x

j

|,

where theminimization isover allvorti esandimage vorti es.

4.4.4 Computation of temperatures

We he k the inverse temperature

β

and reservoir varian e

γ

numeri ally

as-suming ergodi ity. For some fun tion

a(X) :

D → R

d

and an equilibrium

distribution

ρ(X) = exp

h

−β

∗

_H(X)

_˜

i

∇ · ρ(X)a(X) = a(X) · ∇ρ(X) + ρ(X)∇ · a(X)

=

_−β

∗

ρ(X)a(X)

_{· ∇ ˜}

H + ρ(X)

_{∇ · a(X).}

Formallyintegratingoverphasespa e

Z

D

∇ · ρ(X)a(X)dX = − β

∗

Z

D

ρ(X)a(X)

· ∇ ˜

H dX +

Z

D

ρ(X)

∇ · a(X)dX.

(4.28)

The expressionon theleft is zeroifeither

ρ

or

a

is zeroon theboundary

∂

D

ofphasespa e. Theboundaryof

D

onsistsof ongurationsforwhi hatleast

one pointvortex islo ated ontheboundaryof thedis . Su h a onguration

has energy

H

→ −∞

. Likewise, there are points in phase spa e where two

or more point vorti es ollide and the Hamiltonian tendsto

±∞

. TheGibbs

distribution (4.6) anbenormalized onlyfor

β

onanopeninterval[9 ℄:

β

∈

µ −8π

_Γ

2 _M

,

+4π

Γ

2 ¶

(14)

4.4. A thermostated integrator for point vortices

89

To arryouttheintegration(4.28), we hoose

a

fortheform:

a = b/ρ,

ρ(X) = exp [

_{−βH(X)] ,}

where

β

isthedesiredinversetemperatureand

b(X)

issomefun tionwith

b = 0

attheboundaryofthephasespa e. Inthis ase,theexpressionfor

β

∗

simplies

to

0 =

−β

∗

ha · ∇ ˜

H

i + h∇ · ai.

If the ow is ergodi , then the ensemble averages an be repla ed with time

averages

β

∗

=

∇ · a/a · ∇ ˜

H,

andthedisagreementof

β

∗

and

β

servesasasimple he kfornonergodi ity. For

theinnite reservoir,

H = H

˜

, and forthenite reservoir,

H = H + γ

˜

∗

_/β

∗

_H

2

, yielding

0 =

−β

∗

ha · ∇Hi − 2γ

∗

ha · H∇Hi + h∇ · ai.

Choosing two independentfun tions

a

1 = b

1 /ρ

and

a

2 = b

2 /ρ

,where

b

1

and

b

2

areidenti allyzeroon

∂

D

,theseequationsyieldalinearsystemfor

β

∗

and

γ

∗

.

Forourexperimentswe hose

b

1 =

∇H

Y

i

(R

2 _{− |x}

i

|

2 )

|x

i

|

2 ,

b

2 =

∇H

Y

i

(R

2 _{− |x}

i

|

2 )

2 |x

i

|

4 ,

ρ = exp(

_{−βH − γH}

2 ),

where

β

is one of the three inverse temperatures (4.30) and

γ

is either

0

for

theinnitereservoirorthe orrespondingreservoirvarian e(4.31)forthenite

reservoir.

Figure 4.1 illustrates onvergen e of

β

∗

to the values of

β

(4.30) for both

theinniteandnitereservoir,aswellas onvergen e of

γ

∗

to

γ

(4.31)forthe nite reservoir.

1060

1062

1064

1066

1068

1070

−0.01

−0.005

0

0.005

0.01

0.015

0.02 t

β

∗

1060

1062

1064

1066

1068

1070

−0.01

−0.005

0

0.005

0.01

0.015

0.02 t

β

∗

1060

1062

1064

1066

1068

1070

−1

0

1

2

3

4 x 10

−5

t

γ

∗

A

B

C

A

B

C

A

B

C

Figure

_{4.1: Convergence of inverse temperature β and reservoir variance γ for}

high (A), medium (B) and low (C) temperature states. Left: infinite reservoir

β

∗

_{. Middle: finite reservoir β}

∗

_{. Right: finite reservoir γ}

∗

_.

(15)

4.5 Numerical experiments

Forallofthenumeri alexperimentsusingfourstrongvorti es,theinitial

on-guration onsists of point vorti es with ir ulations and positions given by

[9 ℄:

Γ

1 = Γ

3 = 10π,

Γ

2 = Γ

4 =

−10π,

x

1 = (3, 0),

x

3 = (

−3, 0), x

2 = (0, 3),

x

4 = (0,

−3).

Forboththeniteandinnitereservoirthermostatwe hoosenegative,neutral

andpositiveinversetemperatures

β =

{−0.006, −0.00055, 0.01}.

(4.30)

Our hoi e of

β

is su h that itis loseto the theoreti alupper limit in (4.29)

forhigh temperature,and itis loseto thetheoreti allower limitin (4.29)for

lowtemperature.

Thesizeofthereservoir isdened by

γ

. Inthe aseofaninnitereservoir

γ

_{≡ 0}

,fora nitereservoir

γ = β/(

_−2E

0 )

with

E

0 =

{628, 221, −197}.

(4.31)

Inallexperiments,wetake

α = 0.5

and

λ =

√

0.4

.

Weintegrated the thermostated dynami s over the interval

t

∈ [0, T ]

with

T = 12000

using the time transformation (4.27) and xed transformed time

steps

∆s = ∆t

0 /e(X(0))

with

∆t

0 = 0.001

. Thesamplingwasperformedover

the time interval

[T

0 , T ]

with

T

0 = 1500

, to allow de orrelation of the initial

onditions. Theresultingtime series was sampled uniformly in time in y les

of

δt = 0.01

, toprodu ethehistogramsshown inFigures4.64.9.

4.5.1 Ergodicity tests

The extendedmeasure (4.12) isGaussian in the thermostat variable

ζ

. If the

time dynami sisergodi withrespe tto(4.12), weexpe t thetimeseries

ζ(t)

tobenormallydistributed,i.e.

ζ

∈ N (0, α

−1

₎

. Ahistogramofthevaluesof

ζ

is

showninFigure4.2fortheneutral ase

β =

−0.00055

. Thenormaldistribution

ρ(ζ) =

p

α

2π

exp(

−

α

2 ζ

2 )

is also plotted in the gure. The agreement is good,

indi atingergodi itywithrespe tto

ζ

.

As a se ond indi ationof ergodi ity, weplot themotion ofa single vortex

x

1 (t)

in Figure 4.3. The motion appears well-mixed. The density of points

along the traje tory is greater where either the lo al velo ity

˙x

1

or the lo al

timestep

∆t

n

issmall.

4.5.2 Momentum conservation

Thefun tion

s

1 (X)

in(4.23)is hosentopreserveangularmomentum(4.21)of

the strongvortex set underthe thermostated dynami s. Figure 4.4 showsthe

angular momentum

M

ˆ

as a fun tion of time for the three temperatures. We

(16)

4.5. Numerical experiments

91 −6

−4

−2

0

2

4

6

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4 ρ

(ζ

)

ζ

Figure

_{4.2: Distribution of thermostat (dash line), Gaussian fit (solid line).}

−5

0

5 −5

−4

−3

−2

−1

0

1

2

3

4

5 Figure

_{4.3: Motion of a single vortex x}

₁

_{(t) on the interval t}

_{∈ [0, 1000] for}

(17)

0 2000

4000

6000

8000

10000

12000

14000

16000

18000

−7

−6

−5

−4

−3

−2

−1

0

1

2

3 x 10

−10

t

ˆ M

(t

)

A

B

C

Figure

_{4.4: Momentum for high (A), medium (B) and low (C) temperature}

states for finite reservoir size. Infinite reservoir gives a similar behavior.

4.5.3 Temperature effects

Inthis se tionweattempt to reprodu etheexperimentsof Bühlerusing

ther-mostatedlargepointvorti es. We ondu t experimentsusingboththeinnite

reservoir anoni al distribution ((4.22) with

γ

≡ 0

) and the nite reservoir

distribution ((4.22)with

γ

6= 0

).

The time evolution of thekineti energy of strong vorti es is displayed in

Figure 4.5 for both the innite and nite reservoir models, showing that the

thermostatdrivestheenergyevolutiontowardsthedesiredtemperature. Figure

4.6showstheprobabilitydistributionsofthekineti energyofthevorti es. For

thenitereservoirthermostat,themeansandvarian esaresimilarto thoseof

[9 ℄.

Figure4.7displaysthehistogramofdistan es

|x

i

− x

j

|

betweenlike-signed

vorti es. Biasinfavorofsmallseparationsisevidentatnegativetemperatures,

onsistent with Onsager's predi tions. The distributions are very similar to

those obtained byBühler [9℄. For theinnitereservoir model, there isa large

peak in the distribution at

|x

i

− x

j

| ≈ 1

whi h is in onsistent with Bühler's

simulations. Thiso ursbe ausetoomu henergyisdrawnfromthereservoir.

The omparisonisre overedin thenitereservoirmodel.

Figure 4.8 shows the histograms of the distan e between opposite-signed

vorti es. Inthis ase,thereisasomewhatmilderbiastowards loseapproa hes

at negative temperatures, in keeping with Onsager's ideas. The bias is less

(18)

oppo-4.5. Numerical experiments

93 2000

4000

6000

8000

10000

12000

−400

−300

−200

−100

0

100

200

300

400

500

600 t

H

(t

)

Infinite reservoir

2000

4000

6000

8000

10000

12000

−400

−300

−200

−100

0

100

200

300

400

500

600 t

H

(t

)

Finite reservoir

A

B

C

A

B

C

Figure

_{4.5: Time evolution of energy H(t) for infinite (left) and finite (right)}

reservoirs. Inverse temperatures: β = 0.01 (solid line), β =

_{−0.00055 (line with}

x-mark), β =

_{−0.006 (dash line).}

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

A

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

B

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

C

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

A

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

B

−1000

0 −500

0

500 1000

1

2

3

4

5 ρ

(H

)∗

10

3 H

C

Figure

_{4.6: Distribution of energy for high (A), medium (B) and low (C)}

(19)

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

A

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

B

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

C

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

A

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

B

0

2

4

6

8

10

0

0.5

1

1.5

2

2.5 ρ

(|

x

i

−

x

j

|)

∗

10 |x

i

−x

j

|

C

Figure

_{4.7: Interparticle spacing among same-signed vortices for high (A),}

medium (B) and low (C) temperature states. Top: Infinite reservoir size.

Bot-tom: Finite reservoir size.

Again the histograms are in ex ellent agreement with the simulation data of

[9 ℄, for the nite reservoir simulation. For an innite reservoir, the positive