Including dynamical fermions in simulations

following), which has previously been defined in eq. 3.19

F_μ^G(x) ≡2[U_μ(x)F_μ^G(x)]T.A.. (3.44) From equation 3.39 we obtain the explicit expression for F_μ^G in the case of the Wilson gauge action

F_μ^G(^x) = ^β 6

∑

ν =μ(^Vμν^(u)(^x)−^Vμν^(d)(^x)), (3.45) whereV_μν^(u) and V_μν^(d) are defined in eq. 3.40. Theforce magnitudeis a real number defined by

||F^G||² ≡2Tr{F^GF^G†} ≥0, (3.46) and it will be recalled quite often in future considerations, also for different fermion forces that we will consider. ||^F||²is given by the trace of a sum of closed paths, therefore this expression is gauge invariant. As noted in [59], if we for example consider the gauge force in directionμ =0, we have

||F₀^G(x)||²= 2β

3 ₂

∑

Tr{E_iE_j}, (3.47)

where E_i(x) ≡ ¹₂U₀(x)(V_0i^(u)(x)−V_0i^(d)(x))

T.A. is the discretized electric field in lattice units.

3.7 Including dynamical fermions in simulations

So far we have considered a way to simulate the theory that takes into account only the gauge fields dynamics. This corresponds to neglecting the vacuum loops of quarks that are present in the full theory. The QCD partition function which, in addition to the gauge actionS_G[U], also includes fermionic dynamical degrees of freedom is given by the following

Z = [^DU][Dψ][Dψ]^{e−SG(U)−SF(U,ψ,ψ)}. (3.48)

If we consider the case with two degenerate flavors of dynamical quarks, the corresponding fermion action reads

S_F(U,ψ,ψ) =

∑

f=1

ψ_fD(^U)ψ_f, (3.49)

whereD(U) = D₀(U) +mrepresents the massive Dirac operator.

The eigenvalues of the lattice Dirac operator D(^U)are complex and their real parts are not always positive definite. Therefore, if one would like to simulate just a single fermion flavor with the simple approach described abogve, the exponen-tial in eq. 3.48 would not necessarily have to real and positive and in general cannot be interpreted as a conventional probability measure. This is the reason we were from the start considering two degenerate dynamical quark flavors - in this case the Dirac operator enters quadratically in the effective action and the positivity is guaranteed.

Since Grassmann fields ψ_f, ψ_f are very difficult to incorporate directly into a computational simulation, we try to integrate over them in order to obtain the partition function that depends only on the gauge fields and on which the previ-ously discussed algorithms could hopefully be applied. After taking into account the transformation properties of Grassmann variables, we can indeed integrate out the fermion field variables and the partition function now reads

Z= [^DU][^Dψ][^Dψ]^e−SG(U)(detD(^U))². (3.50) For a detailed derivation, see for example [29]. Due toγ5-hermicity of the Dirac operators considered in this work (γ5D = D^†γ5) we have detD = detD^† and in the future we will take advantage of this property by writing det[D²(U)] as det[^D(^U)^D†(^U)].

It is evident that the methods used to simulate pure gauge theory need to be modified to include the effect incorporated in the fermion determinant detD(^U). This determinant is a non-local object, therefore its computation after each change of the gauge field would be numerically extremely expensive. The way out of this problem is to simulate the determinant by introducing thepseudofermion field φ.

This is a bosonic field that has the same number of degrees of freedom as the cor-responding fermionic variable (color, Dirac and space-time). Its introduction was motivated by the analogy between the fermionic and bosonic Gaussian integrals which allows us to express the fermion determinant as a Gaussian integral of the

3.7 Including dynamical fermions in simulations

pseudofermion fields

det[D^†D] = [Dψ][Dψ]e^{−∑2f=1ψfD(U)ψf} ∝ [Dφ][Dφ^†]e^{−φ†(D†(U)D(U))−1φ}, (3.51) where the identity holds up to an irrelevant constant. The partition function with the new bosonic integration variables, for the considered system of two degener-ate fermion flavors, is now proportional to

Z = [DU]e^−SGdet(D^†D) = [DU][Dφ][Dφ]e^−Seff (3.52) S^eff(U,φ,φ^†) =^SG(^U) +^SPF(U,φ,φ^†), (3.53) S_PF(U,φ,φ^†) =−φ^†(D(U)D^†(U))⁻¹φ. (3.54) If the determinant detD(^U) is itself positive, then one may even be able to simulate a single fermion flavor from eq. 3.51. In this case, one approximates an inverseD⁻¹(U)by some operatorTT^†during the molecular dynamics part of the HMC and corrects for this approximation in the accept-reject step at the end of the trajectory. For the approximation of TT^† different polynomial or rational functions can be used. This approach is widely used in simulations of an odd number of flavors or in simulations with degenerate quark masses taken into account dynamically, for example the so-called N_f = 2+1 simulations, which include 2 degenerate massless u,d quarks and a heavier s quark. Since in this work we are interested exclusively in simulating an even number of dynamical fermion flavors with degenerate masses (N_f =2 andN_f =4), we will not discuss this approach further.

The argumentation similar to the one for two flavors holds for applying an arbitraryv number of pseudofermions. For completion, we define the partition function for the theory of N_f fermion flavors, which can be interpreted as a prob-ability measure ifN_fis even

Z = [DU][Dφ][Dφ]e^−SGdet(D^†D)^Nf/2 = [DU][Dφ][Dφ]e^−Seff (3.55)

S^eff =S_G−S_PF, (3.56)

S_PF =−φ₁^†(^DD†)⁻¹φ₁−φ₂^†(^DD†)⁻¹φ2+. . .φ^†_Nf/2(^DD†)⁻¹φ_Nf/2. (3.57) In the above formulae, we have dropped the dependence of the Dirac operator on the gauge field and from now on we will write it only when necessary.

The HMC algorithm may now be applied for simulating the theory with an even number of fermion flavors N_f. The creation of pseudo fermion fields is sim-ple: they can, for example, be created by generating a complex vector χ with

Gaussian distributione^{−χ χ} and then determiningφ = ^Dχ. They do not get up-dated during the HMC trajectory. The modification needed for the procedure discussed in section 3.4 is that now the momenta are generated with probability distribution

P(π) ∝e^12(π,π). (3.58) The HMC steps for updating the gauge fields and momenta follow as for the pure gauge case, but instead of the one in eq. 3.36 the Hamiltonian used here is extended with the pseudofermion action

H(π,U) = ¹

2(π,π) +S_G(U) +S_PF(U,φ,φ^†). (3.59) The force that drives the molecular dynamics evolution now has a part com-ing from the gauge action (discussed in section 3.6) and a part comcom-ing from the pseudofermion action which will be discussed in detail in Chapter 4.

Introducing dynamical fermions into simulations was an important and neces-sary step in lattice QCD simulations for turning them into a precision tool. The effect of considering the fermions dynamically is not negligible and we illustrate it (cf. Figure 3.1) by quoting the figure from [60], which compares the values of several observables computed in pure gauge theory and in the theory with N_f =2+1 dynamical fermions in theasqtadregularization.