Critical behavior for the model of random spatial permutations

Dissertation defense

John Kerl

Department of Mathematics, University of Arizona

March 22, 2010

Outline

• The probability model

• Order parameters and criticality • Markov chain Monte Carlo methods • Finite-size scaling

• The worm algorithm

• Other work

The probability model The probability model: definitions

State space: Ω_{Λ,N} = Λ^{N} × S_{N}, where Λ = [0, L]^{3} with periodic boundary conditions. Point positions: X = (x_{1}, …, x_{N}) for x_{1}, …, x_{N} ∈ Λ. Distance function (short-jump regime with periodic boundary conditions):

∥x − y∥_{Λ} = ∥d∥_{Λ} = min_{n∈ℤ^{3}} ∥d + nL∥.

Hamiltonian, where T = 1/β and r_{ℓ}(π) is the number of ℓ-cycles in π:

H(X, π) = \frac{T}{4} ∑_{i=1}^{N} ∥x_{i} − x_{π(i)}∥_{Λ}^{2} + ∑_{ℓ=1}^{N} α_{ℓ}r_{ℓ}(π).

• The first term discourages long permutation jumps, moreso for higher T.

• The temperature scale factor T/4, not β/4, is surprising but correct for the Bose-gas derivation of the Hamiltonian.

• The second term discourages cycles of length ℓ, moreso for higher α_{ℓ}. These

interactions are not between points, but rather between permutation jumps.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  4 / 28

The probability model: definitions

Fixed point positions (quenched model — includes all simulations done up to the present on the cubic unit lattice with N = L^{3}):

P_{X}(π) = \frac{1}{Y (Λ, X)} e^{−H(X,π)}, Y (Λ, X) = ∑_{σ∈S_{N}} e^{−H(X,σ)}.

Varying positions (annealed model — many theoretical results are available):

P (π) = \frac{1}{Z(Λ, N)} e^{−H(X,π)}, Z(Λ, N) = \frac{1}{N!} ∫_{Λ^{N}} Y (Λ, X) dX.

In either case, we write the expectation of an RV S(π) as 𝔼[S] = ∑_{π∈S_{N}} P (π)S(π).

Feynman (1953) studied long cycles in the context of Bose-Einstein condensation for interacting systems. See also Sütő (1993, 2002), and papers of Betz and Ueltschi.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  5 / 28

The probability model: intuition

What does a typical random spatial permutation actually look like? (Recall H(X, π) = \frac{T}{4} ∑_{i=1}^{N} ∥x_{i} − x_{π(i)}∥_{Λ}^{2} + ∑_{ℓ=1}^{N} α_{ℓ}r_{ℓ}(π).)

• As T → ∞, the probability measure becomes supported only on the identity permutation. Large but finite T: there are tiny islands of 2-cycles, 3-cycles, etc.

• As T → 0, length-dependent terms go to zero. The probability measure approaches

the uniform distribution on S_{N}: all π’s are equally likely. For intermediate T, things get more interesting:

• Lengths of each jump, ∥π(x) − x∥_{Λ}, remain small: empirically, < 3.

• Above a critical temperature T_{c}, all cycles are short: 2-cycles, 3-cycles, etc. T_{c} ≈ 6.87, and positive α terms increase T_{c}.

• Phase transition at T_{c}: below T_{c}, jump lengths remain short but long cycles form.

Order-parameter RVs 1/ξ, f_{I}, f_{M}, f_{W}, f_{S} (below) quantify this. High T, medium but subcritical T, and low T:

Order parameters and criticality Order parameters: 1/ξ, f_{S}, f_{W}, f_{I}, f_{M}

The spatial cycle length and correlation length are

s_{x}(π) = ∑_{j=1}^{ℓ_{x}(π)} ∥π^{j}(x) − π^{j−1}(x)∥_{Λ}  and  ξ = \overline{s}(π) = \frac{1}{N} ∑_{x∈Λ} s_{x}(π).

The winding number of π counts the integer number of wraps of π’s cycles around the 3-torus in each of the three directions:

W(π) = (W_{x}(π), W_{y}(π), W_{z}(π)) = \frac{1}{L} ∑_{i=1}^{N} ∥π(x_{i}) − x_{i}∥_{Λ}
W^{2}(π) = W(π) · W(π) = W_{x}(π)^{2} + W_{y}(π)^{2} + W_{z}(π)^{2}.

The scaled winding number, f_{S}, arises in physics:

f_{S} = \frac{𝔼[W^{2}]T L^{2}}{3N} = \frac{𝔼[W^{2}]T}{3L} .

The fraction of sites in cycles which wind, f_{W}: self-explanatory.

The fraction of sites in long cycles, f_{I}: defined in the dissertation. Intuition matches the name.

The scaled mean longest cycle length: f_{M} = 𝔼[ℓ_{max}]/N.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  8 / 28

Behavior of order parameters as functions of L, T, and α.

1/ξ is right-sided; the rest are left-sided. All order-parameter plots tend to the right as α increases, i.e. ΔT_{c}(α) = \frac{T_{c}(α)−T_{c}(0)}{T_{c}(0)} is positive for small positive α. Goal: quantify ΔT_{c}(α)’s first-order dependence on small α.

Known results and conjectures

Recall H(X, π) = \frac{T}{4} ∑_{i=1}^{N} ∥x_{i} − x_{π(i)}∥_{Λ}^{2} + ∑_{ℓ=1}^{N} α_{ℓ}r_{ℓ}(π). We have the following models:

• Non-interacting model: α_{ℓ} ≡ 0.

• Two-cycle model: α_{2} = α and other cycle weights are zero.

• Ewens model: α_{ℓ} is constant in ℓ.

• General-cycle model: No restrictions on α_{ℓ}.

Known results for points on the continuum (obtained largely using Fourier methods):

• ΔT_{c}(α) is known (to first order in α) for two-cycle interactions (Betz and Ueltschi, CMP 2008) and small cycle weights (Betz and Ueltschi 2008). (This taps into a long and controversial history in the physics literature: see Baym et al., EJP B 2001, or Seiringer and Ueltschi, PRB 2009, for surveys.) The critical (ρ, T, α) manifold relates ρ_{c} to T_{c}.

ρ_{c}(α) ≈ ∑_{ℓ≥1} e^{−α_{ℓ}} ∫_{ℝ^{3}} e^{−ℓ 4π^{2}β∥k∥^{2}} dk = \frac{1}{(4πβ)^{3/2}} ∑_{ℓ≥1} e^{−α_{ℓ}}ℓ^{−3/2}
ΔT_{c}(α) ≈ cρ^{1/3}α, for small α, with c ≈ 2/3 when ρ = 1.

Markov chain Monte Carlo methods Metropolis sampling

The expectation of a random variable S (e.g. f_{W}, f_{M}, f_{I}, f_{S}, ξ) is

𝔼[S] = ∑_{π∈S_{N}} P (π)S(π).

N! grows intractably in N. Instead, estimate expectations by summing over some number M (10^{5} or 10^{6}) typical permutations. The sample mean is now a random variable with its own variance.

The usual technical issues of Markov chain Monte Carlo (MCMC) methods are known and handled in my simulations and dissertation: thermalization time, proofs of irreducibility, aperiodicity, and detailed balance (below), autocorrelation, batched means, and quantification of variance of sample means.

The fundamental Metropolis step (the analogue of single spin-flips for the Ising model) swaps permutation arrows which end at nearest-neighbor lattice sites. This either splits a common cycle, or merges disjoint cycles:

As usual, the proposed change is accepted with probability min{1, e^{−ΔH}}.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  12 / 28

Correctness for the swap-only (SO) algorithm

Detailed balance, i.e. P (π)M(π, π′) = P (π′)M(π′, π) for all π, π′, is easy to prove using standard Metropolis M(π, π′) ∼ min{1, e^{−ΔH}}. Here we prove irreducibility.

Proposition: Any π′ is reachable from any other π using swaps. Proof. Transpositions generate S_{N}. We construct a sequence of (nearest-neighbor) swaps which results in a non-nearest-neighbor swap. We are given π, x, and z. Choose a nearest-neighbor path y_{0} = x, y_{1}, …, y_{n−1}, y_{n} = z. (See figure.) Swaps: (y_{0}, y_{1}), (y_{0}, y_{2}), … (y_{0}, y_{n}); then (y_{n}, y_{n−1}), (y_{n}, y_{n−2}), … (y_{n}, y_{1}). Then π′(x) = π(z), π′(z) = π(x), and π′(y) = π(y) for all other y. □

PSfrag replacements

y_{0} π(y_{0})

y_{1} π(y_{1})

y_{2} π(y_{2})

y_{3} π(y_{3})

y_{3} y_{2} y_{1} y_{0}

π(y_{3}) π(y_{2}) π(y_{1}) π(y_{0})

y_{0} π(y_{0})

y_{1} π(y_{1})

y_{2} π(y_{2})

y_{3} π(y_{3})

y_{3} y_{2} y_{1} y_{0}

π(y_{3}) π(y_{2}) π(y_{1}) π(y_{0})

y_{0} π(y_{0})

y_{1} π(y_{1})

y_{2} π(y_{2})

y_{3} π(y_{3})

y_{3} y_{2} y_{1} y_{0}

π(y_{3}) π(y_{2}) π(y_{1}) π(y_{0})

(5) (4) (3) (2) (1) (6) y_{1} y_{2} z x

Conclusion: given irreducibility, aperiodicity (also easy), and detailed balance, the Gibbs distribution is the invariant (and thus limiting) distribution for the SO chain.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  13 / 28

Conservation of winding number (with probability near 1)

Proposition: If jump lengths are less than L/2, swaps conserve winding number. Proof. Swaps are done on pairs of arrows which end at nearest-neighbor sites. Due to short jump lengths, all four sites in a swap are in the same Euclidean chart lifted off the 3-torus. Thus

W′ − W = \frac{1}{L} ∑_{i=1}^{N} d̃(π′(x_{i}), π(x_{i})) = \frac{1}{L} [d̃(π′(x), π(x)) + d̃(π′(y), π(y))]

= \frac{1}{L} [π′(x) − π(x) + π′(y) − π(y)] = \frac{1}{L} [π(y) − π(x) + π(x) − π(y)] = 0.

□

PSfrag replacements x_{1} x_{2} x_{3} x_{4} x_{5} x_{6} x_{7} x_{8} x_{1} x_{2} x_{3} x_{4} x_{5} x_{6} x_{7} x_{8}

Common Euclidean chart

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  14 / 28

Partial solution: the swap-and-reverse (SAR) algorithm

Figure part 1: a long cycle on the torus almost meets itself in the x direction.

Part 2: after a swap-only step (above), one cycle winds by +1, and the other by −1. Metropolis steps in the short-jump-length regime create winding cycles only in opposite-direction pairs; total W_{x}(π) is still zero.

Part 3: if we reverse one cycle (zero-energy move), W_{x}(π) is now 2. This swap-and-reverse algorithm permits winding numbers of even parity in each of the three axes: one sweep proposes swaps at each lattice site. A second sweep reverses arrows on each cycle in the permutation with probability 1/2.

Using SAR, it still takes a jump length ≈ L/2 — which happens effectively never — to create an odd winding number. Band updates (see dissertation) are one idea; the worm algorithm (below) is another.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  15 / 28

Finite-size scaling Computational results: finite-size scaling method

Raw MCMC data yield S(L, T, α) plots as above, for each order parameter S. Finite-size scaling (see Pelissetto and Vicari, arXiv:cond-mat/0012164, for a survey) determines the critical temperature T_{c}(α).

Define reduced temperature t = \frac{T −T_{c}(α)}{T_{c}(α)} , and correlation length ξ as above.

Hypotheses: (1) At infinite volume, S ∼ | − t|^{ρ} and ξ ∼ |t|^{−ν} (power-law behavior).

(2) Finite-volume corrections enter only through a universal function Q_{S} of the ratio L/ξ:

S(L, T, α) = L^{−ρ/ν}Q_{S}((L/ξ)^{1/ν}) = L^{−ρ/ν}Q_{S}(L^{1/ν}t)

Method:

• Estimate critical exponents ρ, ν via power-law regression on MCMC data plots.

• Plot L^{ρ̂/ν̂}S(L, T, α) as function of T. Since t = 0 at T_{c}(α), these plots for different L cross (up to sampling variability) at T_{c}(α).

• Having estimated ρ̂, ν̂, and T̂_{c}(α), plot L^{ρ̂/ν̂}S(L, T, α) as function of L^{1/ν̂}t̂. This causes all curves to collapse, confirming the FSS hypothesis.

• Regress ΔT̂_{c}(α) on α to estimate the constant c.

Computational results: power-law fit and crossing plots

Raw data vs. power-law fit for 1/ξ and f_{S} with α = 0:

Plots for 1/ξ, f_{I}, and f_{M} show crossing; plots for f_{S} and f_{W} do not. This is most clear at L = 30, 40, 50 where I did M = 10^{6} MCMC samples for T near T_{c}, and most clear in the power-law-fit point of view:

Collapse plots and ΔT̂_{c}(α) ≈ ĉα, given ρ̂’s, ν̂, and T̂_{c}(α)

Omit α = 0.004 since T̂_{c}(α) begins to curve. Omit f_{S} and f_{W} due to non-crossing. Regressing on the (α, ΔT̂_{c}(α)) data, we find T̂_{c}(0) ≈ 6.873 ± 0.006 and ĉ ≈ 0.618 ± 0.086 (2 σ error bars) for Ewens weights on the lattice. For small cycle weights on the continuum, Betz and Ueltschi have T_{c}(0) ≈ 6.625 and c ≈ 0.667. Conclusions: (1) Lattice structure modifies the critical temperature; (2) the α-dependent shift in critical temperature is unaffected.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  19 / 28

The worm algorithm The worm algorithm: intuition

Random-cycle model with p.b.c. has multiple energy minima, indexed by winding numbers. Draw from path-integral Monte Carlo (PIMC) methods in physics. Tunnel through the energy barriers by opening a cycle, modifying it with swap-type moves at its tips, and closing it. Central point: this open cycle, or worm, can wander around the 3-torus (too) freely.

Not permutations anymore? In the figure: nothing ↦ x_{1}, x_{1} ↦ x_{2}, x_{2} ↦ x_{3}, x_{3} ↦ nothing. Rename the nothing to something, called the wormhole point, or w. It has no spatial coordinates and zero distance from any point. We now have π ∈ S_{N+1}, an extended lattice, and an extended random-cycle model.

Same recipe applies as before: (extended) energy function and Metropolis moves; prove correctness. Invent any convenient extended energy for open π’s agreeing with the original energy H for closed π’s (proved next). Sample RVs only on closed π’s.

The worm algorithm: marginality

Proposition: Let S_{N} ↪ S_{N+1} by π(w) = w. Let H, H′ be energy functions on S_{N} and S_{N+1} such that for all π ∈ S_{N}, H(π) = H′(π). Let P, P′, Z, Z′ be as above. Then for all π ∈ S_{N}, P′(π | π ∈ S_{N}) = P (π). Proof.: Let π ∈ S_{N}. By definition of conditional expectation,

P′(π | π ∈ S_{N}) = \frac{P′(π) 1_{S_{N}} (π)}{P′(S_{N})} .

The numerator is Gibbs P for closed permutations, or 0 for open ones:

P′(π) 1_{S_{N}} (π) = \frac{1}{Z′} e^{−H′(π)} 1_{S_{N}} (π) = \frac{1}{Z′} e^{−H(π)} 1_{S_{N}} (π)

since H and H′ agree on closed π’s. The denominator is total probability of closed permutations:

P′(S_{N}) = \frac{1}{Z′} ∑_{π∈S_{N}} e^{−H′(π)} = \frac{1}{Z′} ∑_{π∈S_{N}} e^{−H(π)}.

Since π ∈ S_{N}, the ratio is

\frac{\frac{1}{Z′} e^{−H(π)} 1_{S_{N}} (π)}{\frac{1}{Z′} ∑_{π∈S_{N}} e^{−H(π)}} = \frac{e^{−H(π)} 1_{S_{N}} (π)}{∑_{π∈S_{N}} e^{−H(π)}} = \frac{e^{−H(π)} 1_{S_{N}} (π)}{Z} = P (π).
□.

J. Kerl (Arizona)  Critical behavior for random spatial permutations  March 22, 2010  22 / 28

The worm algorithm: Metropolis moves

PSfrag replacements

x π(x)

w x

π(x)

w

Open at x w.p.

∼ 1 ∧ e^{−ΔH}

π^{−1}(w) π(w)

w π^{−1}(w)

π(w)

w Close w.p.

∼ 1 ∧ e^{−ΔH}

x π(x)

Head swap at x w.p.

∼ 1 ∧ e^{−ΔH} Tail swap at x w.p.

∼ 1 ∧ e^{−ΔH}

x π(x)

π^{−1}(w) w x

π(x) π^{−1}(w)

w x

π(x)

w π(w)

w π(w)

The worm algorithm: fibration and correctness

Key to proving correctness: fibration of S_{N+1} over S_{N}. • N open permutations lie over each closed permutation; fibers partition S_{N+1}.

• Opens and closes stay within fibers.

• Head swaps and tail swaps move across fibers, and furthermore are transitive on fibers.

• Any SO swap can be constructed by an open, a head swap, and a close. (Hence irreducibility via SO, opens, and closes. Aperiodicity and detailed balance: also easy.)

The worm algorithm: stopping time

Good news: examination of random-variable plots for L = 10, comparing SAR to worm, shows that similar results are produced — other than, of course, the winding-number histogram itself.

Problem: The the open worm tips wander around randomly within the L box, and fail to reconnect as L increases. Specifically, histograms show that the distribution of the wormspan ∥π(w) − π^{−1}(w)∥ peaks around L/2.

SAR and worm CPU times are both ∼ aN + bN^{2}. (Shown: L = 5 to 12.) SAR’s b is tiny; worm’s b is not. Interesting L (40-80 or so) are out of reach for the worm algorithm.

Other work Other work

Dissertation items not presented today:

• Precise exposition of the theory of autocorrelation estimators for exponentially correlated Markov processes. Precise quantification of the advantages and non-advantages of batched means.

• Mean length of longest cycle as a fraction of the number of sites in long cycles recovers work of Shepp and Lloyd (1966) for non-spatial uniform permutations.

For the future:

• Use varying (annealed) point positions on the continuum. This samples from the true point distribution.

• Replace cycle-weight interactions in the Hamiltonian with those derived from the true Bose-gas model. Analytical as well as simulational work is needed in order to make this computationally tractable.

Thank you for attending!