The HMC algorithm - Non-Hermitian polynomial hybrid Monte Carlo

The basic idea of HMC is to use the molecular dynamics evolution as transi-tion functransi-tion to evolve from one configuratransi-tion to the next one being accepted or rejected according to a generalized Metropolis step (cf. Figure 5.1).[45] For QCD a configuration is given by the gauge-fieldUµ(x) consisting of aSU(3)c

matrix for each Dirac index µat each lattice site x.

integrating EoM Molecular dynamics

Initial steps

•Generate momentaΠµ

•Generate pseudo-fermionsφ

•Compute HamiltonianHold

Final steps

•Compute HamiltonianHnew

•Accept new configuration ifη < exp{Hnew−Hold} otherwise keep old configuration

Figure 5.1. Schematic flow of the HMC algorithm.

Like in the MD approach we introduce again a fictitious Monte Carlo time τ and guide the evolution by canonical equations. To build the Hamiltonian we generate Gaussian momenta Πµ(x) as conjugate variables toUµ(x). Hence the Π_µare traceless and Hermitian elements of theSU(3)_cgauge group. The

The HMC algorithm 55

Hamiltonian is given by H= 1

2Π²_µ(x) +Seff(Uµ, φ), (5.5)

and evolved according to the equations of motion (EoM) φ˙µ= δ{H}

δ{Πµ} = Πµ; ˙Πµ=−δ{H}

δ{Uµ} =−δ{S_eff}

δ{Uµ}. (5.6) The effective action Seff in case of even-odd preconditioning is the sum of (1.4), (3.17) and (3.18)

Seff(Uµ, φ) =SG(Uµ) +Sdet(Uµ) +Sb(Uµ, φ). (5.7) of the contributions from the gauge-field (G), the factorized part of the de-terminant (det) and the bosonic part (b) which in addition depends on the pseudo-fermion field φ.

The important difference to the MD approach comes at the end of a trajectory when the accept/reject-step is performed. Therefore we compute the energy-difference ∆H between the new Hamiltonian at the end and the old Hamiltonian at the beginning of the trajectory. The probability to accept a new configuration is given by

P_A((U_µ,Π_µ)→(U_µ⁰,Π⁰_µ)) = min{1,exp{−∆H}}, (5.8) with ∆H=Hnew−Hold.

The constructed Hamiltonian serves two distinct tasks: on the one hand it enters in the accept/reject step (acceptance Hamiltonian), on the other hand it guides the molecular dynamics evolution (guidance Hamiltonian).[45] The acceptance Hamiltonian defines the equilibrium distribution to be simulated.

Due to the accept/reject step the algorithm becomes exact accounting for possible rounding or discretization errors occurring e.g. by the numerical in-tegration of the equation of motion. In principal the guidance Hamiltonian entering in the EoM can differ from the acceptance Hamiltonian allowing scope for optimization. Only if both Hamiltonian agree we find energy con-servation (∆H = 0) in the limit of perfect integration.

To let the system evolve according to the EoM we have to compute the variation of the effective action (5.7) in terms of an infinitesimal change of the gauge link δ{Uµ(x)}

δ{Seff}=^X

x,µTr ⁿFµ(x)δ{Uµ(x)}+F_µ^†(x)δ{Uµ(x)^†}^o. (5.9)

The quantity Fµ(x) can be considered as “force” arising for a gauge link varied. It can be split into the contributions

Fµ(x) = Vµ(x) +F_µ^det(x) +F_µ^b(x), (5.10) form the pure gauge part (V_µ), the determinant and the bosonic contribution.

Varying the contribution of the gauge field S_G to the action we obtain the gauge force

V_µ =−β 6

Tr δ{U_P +U_P^†}. (5.11) Next we focus on the determinant force which depends on our choice of preconditioning. Assuming both Mee and Moo to be factorized we are discussing the case of symmetric even-odd preconditioning, where S_det is given by

Sdet=−2 [ln det{Mee}+ ln det{Moo}]. (5.12) Calculating the variation of (5.12) with respect to a link U_µ we find using det{A}= exp{Tr lnA}

δ{S_det}=−2δ{Tr ln(M_ee) + Tr ln(M_oo)}

=−2 Tr^hM^-1_eeδ{Mee}+ M^-1_ooδ{Moo}ⁱ. (5.13) Let us now look at one particular link (pointing from x to x +µ) to be varied and consider like Jansen and Liu all “clover leaves” which contain this link.[32] We find the following four diagrams forµ6=ν as shown in Fig. 5.2.

x x+µ x x+µ x x+µ x x+µ

x+ν x+ν

x−ν x−ν

xodd

xodd xeven xeven

Figure 5.2. Diagrams contributing toF_µ^det(x) with the triangles denoting “even” inser-tions.

The HMC algorithm 57 To compute the contribution by varying the link Uµ(x) we have to com-pute the staples starting atx+µ and follow the arrows. At positions of the black triangles a (3×3)-color matrix has to be inserted which for an even insertion point is given by

(N)x = TrDirac

hiσµνM^-1_ee(x)ⁱ, (5.14) while for an odd point M^-1_ee is replaced by M^-1_oo. Finally, we have to sum over all directionsµ6=ν and multiply the sum by −κc_sw/4 i.e.

F_µ^det(x) = −κcsw

µ6=ν

(diagrams in Fig.5.2). (5.15) The inversion of Mee and Moo is needed locally for all lattice sites x and computed e.g. exactly by applying the Householder triangularization.[36]

The last contribution, named bosonic force, comes from the variation of the pseudo-fermion action

Sb =φ^†QˆQˆ^†⁻¹φ, (5.16) where the pseudo-fermion fields φ are generated from a Gaussian random vector η by multiplying Q to achieve a distribution according to (5.16)

φ = ˆQη. (5.17)

Due to even-odd preconditioning ˆQ is a mapping from odd sites on odd sites with the even sites only created/used in between. Hence η and φ are vectors on odd sites only. Computing the variation of (5.16) we find applying δ{A⁻¹}=−A⁻¹δ{A}A⁻¹

δ{Sb}=−φ^†Qˆ^†−¹^hδ{Qˆ^†}Qˆ^†−¹+ ˆQ⁻¹δ{Q}ˆ ⁱQˆ⁻¹φ. (5.18) The variation of the symmetric preconditioned ˆQ is given by

δ{Qˆ^S}=γ5δ{M^-1_ooMoeM^-1_eeMeo}

=−γ₅M^-1_oo[^M^oe^M^-1ee,1I]^h_δ_{_M⁰_oe_} ^δ^{^M₀^eo^}^{i h}^M^-1^ee₁_I^M^eoⁱ

+γ5M^-1_oo[^M^oe^M^-1ee,1I]^h^δ^{^M₀^ee^}_δ_{_M⁰_oo_}^{i h}_M^-1_oo¹^I_M_oeⁱM^-1_eeMeo (5.19) δ{Qˆ^S^†}=γ5δ{MoeM^-1_eeMeoM^-1_oo}

=−γ₅[^M^oe^M^-1ee,1I]^h_δ_{_M⁰_oe_} ^δ^{^M₀^eo^}^{i h}^M^-1^ee₁_I^M^eoⁱM^-1_oo

+γ₅M_oeM^-1_ee[¹^I^,^M^eo^M^-1oo]^h^δ^{^M₀^ee^}_δ_{_M⁰_oo_}^{i h}^M^-1^ee₁_I^M^eoⁱM^-1_oo, (5.20)

x x+µ x x+µ x x+µ x x+µ x+ν

x+ν

x+ν x+ν

x−ν x−ν

Figure 5.3. The contributing diagrams to F_µ^sw(x).

allowing us to write (5.18) as

δ{S_b}=−X^†γ₅^h_δ_{_M⁰_oe_}^δ^{^M₀^eo^}ⁱY +X^†γ₅^h^δ^{^M₀^ee^} _δ_{_M⁰_oo_}ⁱZ+ H.c. (5.21) with

X =^h⁻^M^-1^ee₁_I^M^eoⁱM^-1_ooQˆ^SQˆ^S^†⁻¹φ (5.22) Y =^h⁻^M^-1^ee₁_I^M^eoⁱQˆ^S⁻¹φ (5.23) Z =^h_M-1⁻¹^I

ooMoe

iM^-1_eeM_eoQˆ^S⁻¹φ. (5.24)

The force F_µ^b splits into one contribution from the hopping terms and one from the Sheikholeslami-Wohlert terms. The hopping force for an even x is given by¹

F_µ^hop(x) =κTrDirac

nγ5(1−γµ)^hYo(x+µ)⊗X_e^†(x)

+Xo(x+µ)⊗Y_e^†(x)^io. (5.25) For the force F_µ^sw resulting from the Sheikholeslami-Wohlert terms we have to consider like for the determinant contribution the staples containing the link pointing form xtox+µas they are shown in Fig. 5.3 and insert for the black boxes the matrix in color space [32]

()x = TrDirac

nγ5σµνZ(x)⊗X(x)^†+ H.c.^o. (5.26) The force F_µ^sw is then given by summing over all diagrams for µ 6= ν and multiplying the result withiκcsw/8.

1Forxoddeis to be replaced by oand vice versa.

Integrator 59

Hence the bosonic force is the sum of both parts F_µ^b(x) =κTrDirac

nγ5(1−γµ)^hYo(x+µ)⊗X_e^†(x) +Xo(x+µ)⊗Y_e^†(x)^io + iκc_sw

µ6=ν

(diagrams in Fig. 5.3). (5.27)

Comparing the three contributions to the total force we find that the computation of the bosonic force is most expensive, but the gauge force has the largest impact. One can account for those differences by improving the integration scheme as is discussed in the next section and demonstrated in Chapter 7.

Im Dokument Non-Hermitian polynomial hybrid Monte Carlo (Seite 60-65)