QCD partition function \[ \mathcal{Z}\ =\ \int\mathcal{D}\bar\psi\,\mathcal{D}\psi\,\mathcal{D}U\, e^{-S[\bar\psi,\psi,U]} \]

• Euclidean lattice action \(S\)

• for definitness we assume \(N_f=2\) mass degenerate Wilson type quark flavours
  • can all be generalised for different \(N_f\)-values

• fermion fields Grassmann valued, difficult to represent on a computer

\(\Rightarrow\quad\) integrate out the fermion fields

\(\Rightarrow\quad\) integrate out the fermion fields

\[ \mathcal{Z}\ =\ \int\,\mathcal{D}U\,[\det(D)]^2\, e^{-S_g[U]} \]

and represent the determinant stochastically

\[ \mathcal{Z}\ =\ \int\mathcal{D}\phi^\dagger\,\mathcal{D}\phi\,\mathcal{D}U\, e^{-\phi^\dagger\frac{1}{D[U]^2}\phi - S_g[U]} \]

• with pure gauge action \(S_g\)

• \(D\) the Lattice Dirac operator

• \(\phi^\dagger,\phi\) so-called pseudo-fermion (bosonic) fields

Traded Grassmann-valued fields with a highly non-local term \[ S_\mathrm{eff} = S_\mathrm{PF} + S_g\,,\qquad S_\mathrm{PF}\ =\ \phi^\dagger\,D[U]^{-2}\,\phi \] due to the matrix inverse!

\(\Rightarrow\quad\) local Metropolis updates unrealistic

• any local change requires the evaluation of \(S_\mathrm{PF}\)

random global updates?

• acceptance rate negligible for non-infinitesimal updates

Recall the steps for the local Metropolis algorithm

• probability density \(g(x), x\in\mathbb{R}^n\) (later \(g\equiv e^{-S}\))

Starting from initial \(x\)

1. generate proposal \(x_i'\) for \(x_i\) with random \(i\)
2. accept/reject \(x'\) with probability \[ P_\mathrm{acc}\ =\ \min\left\{1,\ \frac{g(x')}{g(x)}\right\} \]
3. restart from 1. for another random \(i\)-value

• proposal needs to be ergodic

• Metropolis algorithm fulfils detailed balance \[ g(x)\ P_\mathrm{acc}(x\to x') = g(x')\ P_\mathrm{acc}(x'\to x) \]

\(\Rightarrow\quad\) algorithm generates a Markov chain

• this works well for pure gauge theory
  action local

• again, with fermions this is no longer the case

Goal: smartly replace the local update by a global one!

Introduce an artificial Hamiltonian \[ \mathcal{H}(p, x)\ =\ \frac{1}{2}\sum p^2 - \ln(g(x)) \] with canonical momenta \(p\in\mathbb{R}^n\) and a one-paramter familty of deterministic maps \[ \varphi^\tau(p, x): (p, x) \to (p', x') \] with \(\tau\in\mathbb{R}\) which is reversible, i.e. for all \(p, x\) \[ \varphi^\tau(p, x)=(p', x')\quad \Rightarrow\quad \varphi^\tau(-p',x') = (-p, x) \]

Then the HMC algorithm is defined by the following steps

Starting from initial \(x\in\mathbb{R}^n\)

1. heatbath: \(p\ \sim\ \mathcal{N}_{0,1}\)

2. set \((p', x') = \varphi^\tau(p, x)\) with \(\tau\neq 0\)

3. accept/reject \((p', x')\) with probability \[ P_\mathrm{acc}\ =\ \min\{1, \exp(-\Delta\mathcal{H})\}\,,\quad \Delta\mathcal{H} = \mathcal{H}(p,x) - \mathcal{H}(p', x') \]

4. restart at 1.

Don’t need to know the normalisation: let \(g(x)\propto f(x)\)

The HMC transition probability \[ \begin{split} \mathbb{P}_\mathrm{HMC}(x, x') = \int d p\ &\mathcal{N}_{0,1}(p)\, \delta(x'-\varphi^\tau_x(p, x))\\ &\mathbb{P}_\mathrm{acc}[p, x \to \varphi^\tau(p,x)]\,.\\ \end{split} \] We can introduce an integral over \(p'\) as well, which gives \[ \begin{split} \mathbb{P}_\mathrm{HMC}(x, x') = \int d p\,d p'\ &\mathcal{N}_{0,1}(p)\, \delta((p',x')-\varphi^\tau(p, x))\\ &\mathbb{P}_\mathrm{acc}[p, x \to p',x']\,.\\ \end{split} \]

Due to the Metropolis accept/reject step, we have detailed balance with respect to \(\exp(-\mathcal{H})\), and thus \[ \begin{split} e^{-\mathcal{H}(p, x)}\,&\mathbb{P}_\mathrm{acc}[p, x \to p',x']\\ & = e^{-\mathcal{H}(p', x')}\,\mathbb{P}_\mathrm{acc}[p', x' \to p,x]\,.\\ \end{split} \] Moreover, \[ f(x)\, \mathcal{N}_{0,1}(p) = C e^{-\mathcal{H}(p, x)} \] with a normalisation factor \(C\).

Using these properties and the reversibility of \(\varphi^\tau\) it follows \[ \begin{split} f(x)\, &\mathbb{P}_\mathrm{HMC}(x, x') \\ = C\int &d p\,d p'\, e^{-\mathcal{H}(-p', x')}\, \delta((-p,x)-\varphi^\tau(-p', x'))\\ &\quad\mathbb{P}_\mathrm{acc}[-p', x' \to -p,x]\,,\\ \end{split} \] where we have also used that \(\mathcal{H}\) as well as the acceptance probability are agnostic to the sign of \(p\).
Since \(d pd p' = d(-p)d(-p')\), detailed balance follows.

… so, what about this \(\varphi^\tau\)??

The only option I know is the map \(\varphi^\tau\) induced by Hamilton’s EOM \[ \begin{split} \dot p_i \equiv \frac{d p_i}{d\tau} &= -\frac{\partial\mathcal{H}}{\partial x_i}\,,\qquad \dot x_i \equiv \frac{d x_i}{d\tau} &= \frac{\partial\mathcal{H}}{\partial p_i}\,.\\ \end{split} \]

• parameter \(\tau\) becomes trajectory length
• also called Molecular Dynamics (MD) update

• if Hamilton’s EOM solved exactly, \(\mathcal{H}\) is conserved
  \[ \Rightarrow\quad\Delta\mathcal{H}=0 \] \(\Rightarrow\quad\) 100% acceptance!

• in practice: need to solve EOMs numerically

• use a symplectic, reversible integration scheme
  like for instance the Leap-Frog integration scheme

• those conserve a shadow Hamiltonian \(\mathcal{H}_s\neq\mathcal{H}\)
  as correction in powers of the integration step size \(\Delta\tau\)

Perform a \(\chi^2\) fit using the HMC \[ \mathcal{H} = \frac{1}{2}(p_a^2 + p_b^2) + \chi^2\,,\qquad \chi^2 = \frac{1}{2}\frac{(y - at - b)^2}{\Delta y^2} \]

Derivatives \[ \frac{d\mathcal{H}}{da} = -t\frac{y-at-b}{\Delta y^2}\,,\qquad\frac{d\mathcal{H}}{db}=-\frac{y-at-b}{\Delta y^2} \]