The LHC is making waves

The LHC is making waves on the web already before it is even being switched on. A rather unusual concept of what experimentalists really care about has been featured in this comic, but it is interesting to know that there are people who are genuinely worried that the LHC might "give helicopters cancer" or even swallow up the entire solar system. Peter Steinberg at the US LHC blog has more about them, including a link to this fine example of paranoia in action -- or is it intentional misinformation, or maybe satire? Of course, a well-known thriller writer who likes to claim his fictions as truth has written about CERN as an antimatter factory potentially useful to therrorists looking for an antimatter bomb, so what wonder is it if some members of the woefully uninformed public are willing believe this kind of stuff? More effective outreach is clearly needed.

Word on the arXiv

The arXiv have announced that they now support submissions of "Microsoft Word DOCX or other OOXML (Office Open XML) document[s]". While I am perfectly aware that high-energy physicists (or indeed any kind of physicists) are not the only users of the arXiv, and that usage of TeX is not terribly common outside the physics/mathematics field (though I know a few philosophers and economists, and even one historian, who were won over by the superior look of texts typeset in LaTeX), I find this a little worrying, especially given that the arXiv acknowledges support from the Microsoft Technical Computing Initiative. What worries me is the possibility that this might be the first step towards a less open information architecture at the arXiv, and by implication in the high-energy physics communications sector. Will Microsoft try to gain a foothold, leading to the eventual establishment of their "open" (not) formats as the only accepted submission and download format? One sincerely hopes not.

arXiv catchup

I have been too lazybusy recently to blog anything. However, in the spirit of the day, I'd like to share a romantic little poem extolling the nonabelian nature of strong attraction:

Roses are red, violets are blue
quarks come in colours, and so does glue.

No, I won't give up physics and become a card designer for Hallmark, don't worry. But after softening your hearts with this touching verse, I'd like to blog about some rather old stuff, which I hope hasn't gone stale in the meantime.

One paper on the arXiv that struck me as interesting in the last couple of months was this paper by Jeffrey Mandula (of Coleman-Mandula No-Go fame), who discusses the consequences of Lüscher's nonlinear realisation of chiral symmetry for Ginsparg-Wilson fermions. We recall that this symmetry can be written in two inequivalent ways by putting the phase factor e^iαγ₅ either on the quark field ψ or its conjugate $$\bar{\psi}$$. The crucial fact that Mandula points out is that both of these are independent symmetries of the lattice theory, and they don't commute! Hence, we have to look for the symmetry algebra generated by them, which turns out to be infinite-dimensional. Hence the lattice symmetry has an infinite number of conserved currents, a structure quite different from the continuum theory. However, it would really appear that the differences between any two of these lattice currents are just lattice artifacts of order a or higher that should disappear in the continuum limit, if the latter is properly defined. So some of the objections that the paper raises are likely a lot less serious than stated (especially the non-locality exhibited for free overlap fermions [eq. (38)] goes away once one realises that the continuum limit must be taken with the negative mass s constant in lattice units), but it appears that Ginsparg-Wilson fermions may have their own set of problems beyond just being expensive to simulate. Any comments on this from Ginsparg-Wilson specialists would be of great interest.

Another interesting paper was this one by Mike Creutz who proposed a new fermion discretisation based on features of the electronic structure of graphene. Apparently the low electronic excitations of a grpahene layer are described by the massless Dirac equation, and a lattice model based on this (by reducing the links in one of the three graphene hexagonal directions to points, and rescaling eveything to make the lattice rectangular again) exploits this to achieve the minimum number (two) of doublers permitted in an conventional chiral lattice theory by the Nielsen-Ninomiya theorem, and this construction can be extended to four dimensions and gauged to get a lattice discretisation of QCD with two light quark flavours. This was quickly followed up by a similar proposal for a minimally-doubling quark action, and by this paper which shows that any minimally-doubling chiral lattice theory necessarily has to break either of the discrete symmetries P or T such that their product PT is broken; this allows the generation of additional (relevant) dimension 3 operators that have to be removed by fine-tuning, precluding the use of minimally-doubling chiral actions in practice (unless some additional non-standard symmetry should conspire to do that fine-tuning itself, a possibility hinted at in the conclusion).

Spammers strike

My apologies to our readers for the very nasty spam post that headed this blog for the last few hours. Somehow a spammer managed to break into Matthew's unused Blogger account and used it to post his garbage (in giant bold italics). I have suspended Matthew's account from posting at this blog for the time being, and I hope that this will put an end to the spammers' incursion into this forum. Apologies to all offended parties again.

Algorithms for dynamical fermions -- Hybrid Monte Carlo

In the previous post in this series parallelling our local discussion seminar on this review, we reminded ourselves of some basic ideas of Markov Chain Monte Carlo simulations. In this post, we are going to look at the Hybrid Monte Carlo algorithm.

To simulate lattice theories with dynamical fermions, one wants an exact algorithm that performs global updates, because local updates are not cheap if the action is not local (as is the case with the fermionic determinant), and which can take large steps through configuration space to avoid critical slowing down. An algorithm satisfying these demands is Hybrid Monte Carlo (HMC). HMC is based on the idea of simulating a dynamical system with Hamiltonian H = 1/2 p² + S(q), where one introduces fictitious conjugate momenta p for the original configuration variables q, and treats the action as the potential of the fictitious dynamical system. If one now generates a Markov chain with fixed point distribution e^-H(p,q), then the distribution of q ignoring p (the "marginal distribution") is the desired e^-S(q).

To build such a Markov chain, one alternates two steps: Molecular Dynamics Monte Carlo (MDMC) and momentum refreshment.

MDMC is based on the fact that besides conserving the Hamiltonian, the time evolution of a Hamiltonian system preserves the phase space measure (by Liouville's theorem). So if at the end of a Hamiltonian trajectory of length τ we reverse the momentum, we get a mapping from (p,q) to (-p',q') and vice versa, thus obeying detailed balance: e^-H(p,q) P((-p',q'),(p,q)) = e^-H(p',q') P((p,q),(-p',q')), ensuring the correct fixed-point distribution. Of course, we can't actually exactly integrate Hamilton's in general; instead, we are content with numerical integration with an integrator that preserves the phase space measure exactly (more about which presently), but only approximately conserves the Hamiltonian. We make the algorithm exact nevertheless by adding a Metropolis step that accepts the new configuration with probability e^-δH, where δH is the change in the Hamiltonian under the numerical integration.

The Markov step of MDMC is of course totally degenerate: the transition probability is essentially a δ-distribution, since one can only get to one other configuration from any one configuration, and this relation is reciprocal. So while it does indeed satisfy detailed balance, this Markov step is hopelessly non-egodic.

To make it ergodic without ruining detailed balance, we alternate between MDMC and momentum refreshment, where we redraw the fictitious momenta at random from a Gaussian distribution without regard to their present value or that of the configuration variables q: P((p',q),(p,q))=e^{-1/2 p'2}. Obviously, this step will preserve the desired fixed-point distribution (which is after all simply Gaussian in the momenta). It is also obviously non-ergodic since it never changes the configuration variables q. However, it does allow large changes in the Hamiltonian and breaks the degeneracy of the MDMC step.

While it is generally not possible to prove with any degree of rigour that the combination of MDMC and momentum is ergodic, intuitively and empirically this is indeed the case. What remains to see to make this a practical algorithm is to find numerical integrators that exactly preserve the phase space measure.

This order is fulfilled by symplectic integrators. The basic idea is to consider the time evolution operator exp(τ d/dt) = exp(τ(-∂_qH ∂_p+∂_pH ∂_q)) = exp(τh) as the exponential of a differential operator on phase space. We can then decompose the latter as h = -∂_qH ∂_p+∂_pH ∂_q = P+Q, where P = -∂_qH ∂_p and Q = ∂_pH ∂_q. Since ∂_qH = S'(q) and ∂_pH = p, we can immediately evaluate the action of e^τP and e^τQ on the state (p,q) by applying Taylor's theorem: e^τQ(p,q) = (p,q+τp), and e^τP = (p-τS'(q),q).

Since each of these maps is simply a shear along one direction in phase space, they are clearly area preserving; so are all their powers and mutual products. In order to combine them into a suitable integrator, we need the Baker-Campbell-Hausdorff (BCH) formula.

The BCH formula says that for two elements A,B of an associative algebra, the identity

log(e^Ae^B) = A + (∫₀¹ ((x log x)/(x-1))_{{x=e^{ad A}e^{t ad B}}} dt) (B)

holds, where (ad A )(B) = [A,B], and the exponential and logarithm are defined via their power series (around the identity in the case of the logarithm). Expanding the first few terms, one finds

log(e^Ae^B) = A + B + 1/2 [A,B] + 1/12 [A-B,[A,B]] - 1/24 [B,[A,[A,B]]] + ...

Applying this to a symmetric product, one finds

log(e^{1/2 A}e^Be^{1/2 A}) = A + B + 1/24 [A+2B,[A,B]] + ...

where in both cases the dots denote fifth-order terms.

We can then use this to build symmetric products (we want symmetric products to ensure reversibility) of e^P and e^Q that are equal to e^τh up to some controlled error. The simplest example is

(e^{δτ/2 P}e^{δτ Q}e^{δτ/2 P})^τ/δτ = e^τ(P+Q) + O((δτ)²)

and more complex examples can be found that either reduce the order of the error (although doing so requires one to use negative times steps -δτ as well as positive ones) or minimize the error by splitting the force term P into pieces P_i that each get their own time step δτ_i to account for their different sizes.

Next time we will hear more about how to apply all of this to simulations with dynamical fermions.

Algorithms for dynamical fermions -- preliminaries

It has been a while since we had any posts with proper content on this blog. Lest my readers become convinced that this blog has become a links-only intellectual wasteland, I hereby want to commence a new series on algorithms for dynamical fermions (blogging alongside our discussion seminar at DESY Zeuthen/Humboldt University, where we are reading this review paper; I hope that is not too lazy to lift this blog above the waste level...).

I will assume that readers are familiar with the most basic ideas of Markov Chain Monte Carlo simulations; essentially, one samples the space of states of a system by generating a chain of states using a Markov process (a random process where the transition probability to any other state depends only on the current state, not on any of the prior history of the process). If we call the desired distribution of states Q(x) (which in field theory will be a Boltzmann factor Z^-1e^-S(x)), and the probability that the Markov process takes us to x starting from y P(x,y), we want to require that the Markov process keep Q(x) invariant, i.e. Q(x)=Σ_y P(x,y) Q(y). A sufficient, but not necessary condition for this is the the Markov process satisfy the condition of detailed balance: P(y,x)Q(x)=P(x,y)Q(y).

The simplest algorithm that satisfies detailed balance is the Metropolis algorithm: Chose a candidate x at random and accept it with probability P(x,y) = min(1,Q(x)/Q(y)), or else keep the previous state y as the next state.

Another property that we want our Markov chain to have is that it is ergodic, that is that the probability to go to any state from any other state is non-zero. While in the case of a system with a state space as huge as in the case of a lattice field theory, it may be hard to design an ergodic Markov step, we can achieve ergodicity by chaining several different non-ergodic Markov steps (such as first updating site 1, then site 2, etc.) so as to obtain an overall Markov step that is ergodic. As long as each substep has the right fixed-point distribution Q(x), e.g. by satisfying detailed balance, the overall Markov step will also have Q(x) as its fixed-point distribution, in addition to being ergodic. This justifies generating updates by 'sweeping' through a lattice point by point with local updates.

Unfortunately, successive states of a Markov chain are not really very independent, but in fact have correlations between them. This of course means that one does not get truly independent measurements from evaluating an operator on each of those states. To quantify how correlated successive states are, it is useful to introduce the idea of an autocorrelation time.

It is a theorem (which I won't prove here) that any ergodic Markov process has a fixed-point distribution to which it converges. If we consider P(x,y) as a matrix, this means that it has a unique eigenvalue λ₀=1, and all other eigenvalues λ_i (|λ_i+1|≤|λ_i|) lie in the interior of the unit circle. If we start our process on a state u=Σ_i c_iv_i (where v_i is the eigenvector belonging to λ_i), then P^Nu = Σ_i λ_i^N c_iv_i = c₀v₀ + λ₁^Nc₁v₁ + ..., and hence the leading deviation from the fixed-point distribution decays exponentially with a characteristic time N_exp=-1/log|λ₁| called the exponential autocorrelation time.

Unfortunately, we cannot readily determine the exponential autocorrelation time in any except the very simplest cases, so we have to look for a more accessible measure of autocorrelation. If we measure an observable O on each successive state x_t, we can define the autocorrelation function of O as the t-average of measurements that are d steps apart: C_O(d)=<O(x_t+d)O(x_t)>_t/<O(x_t)²>_t, and the integrated autocorrelation time A_O=Σ_d C_O(d) gives us a measure of how many additional measurements we will need to iron out the effect of autocorrelations.

With these preliminaries out of the way, in the next post we will look at the Hybrid Monte Carlo algorithm.

arXiv API

Via Jacques Distler: The arXiv now has an API intended to allow web application developers access to all of the arXiv data, search and linking facilities. They have a Blog and a Google group about it, as well. Anybody wants to guess when we'll see a "My arXiv papers" application for Facebook?

Warp processors

Via Cosmic Variance, this sounds pretty cool.

Brief post

This is just a brief note saying that I am still alive and still blogging here, but that I have been too busy (moving across the Atlantic, settling into my new job here in Zeuthen, finding a flat in Berlin, and so forth) to write anything since my last post.

To catch up with the aftermath of the lattice meeting: the slides for the plenary and parallel talks are now online (just click on the "pdf" link next to each talk), as are some photos (including one showing me looking into the empty air above, thinking deeply about what to write here). The proceedings are also in progress (the deadline for submissions is 1st October, so expect a flurry of preprints on the lattice arXiv this week).

I'll be back with post that have actual content at some point, but don't expect to much for the next month or so.