Lattice QCD with QCDLAB
Artan Borici

An efficient explicit full-discrete scheme for strong approximation of stochastic Allen-Cahn equationT1
    Xiaojie Wang

[label=e1]x.j.wang7@csu.edu.cn


    December 30, 2023
=======================================================================================================





§ INTRODUCTION


Quantum Chromodynamics (QCD) is the theory of strong interactions. QCD has an ultraviolet fixed point at vanishing coupling constant, a property which was first demonstrated in the perturbative formulation by Gross and Wilczek as well as by Politzer <cit.>. A year later, Wilson was able to formulate QCD non-perturbatively <cit.>. He showed, that in the strong coupling regime, QCD is confining, meaning that the potential between two static charges grows linearly with the separation of charges. Later that year, Kogut and Sussking extended the non-perturbative formulation in the Hamiltonian formalism <cit.>. It was immediately clear that a direct evaluation of QCD path integral was only possible using Monte Carlo simulations. Creutz was the first to show numerically that the weak and strong regimes are in the same phase in four dimensions <cit.>. Since then, lattice QCD has grown into a separate numerical discipline and has delivered results of growing accuracy. This development was possible from the exponential increase of computing power and more efficient algorithms.

In this paper we deal with the basic technology at the bottom of any contemporary lattice computation without going into the details that make lattice QCD confront experiment as well as predict physics beyond the Standard Model. Lattice QCD is a collaborative project, and as such, may not be brought into one review paper without missing a single contribution. Here we profit from the QCDLAB programs which is a small set of short programs that allows one to ilustrate the basic properties of QCD without getting bogged down into the details of advanced computing technology and associated software and algorithms. In contrast to other sotware, QCDLAB maps linear operators of QCD to linear operators of the GNU Octave language <cit.>. Although GNU Octave is an interpreted language, linear operators are precompiled. This property enables very efficient coding as well as minimal run times.

However, GNU Octave is a one-threaded software and runs in one computing core only. Therefore, QCDLAB usage is limited to moderate lattices. It is possible however to include multi-threaded C++ libraries such that the programs run in multiple cores. Writing dedicated libraries of this sort will drive the QCDLAB project out of the original aim of keeping the programing effort small. Nonetheless, Octave is a language in development and is likely to include in the future multi-threaded linear algebra libraries.

In summary, QCDLAB serves three purposes: teaching, learning as well as algorithm prototyping. The latter helps developing a complex software by testing the basic idea of a new algorithm on GNU Octave using QCDLAB codes. QCDLAB programs, version 2.0, as well as this document are available at:





It is licensed under the GNU General Public License v3. The present document serves as a user guide of QCDLAB as well as an illustration of basic calculations in lattice QCD.



§ QCD DATA


The space-time world in lattice QCD is taken to be a four dimesnional regular lattice. At each lattice site go out four directed links, as ilustrated below in the case of two dimensions. The lattice sites are numberd in a lexicographical order.




    < g r a p h i c s >




The basic degree of freedom on the lattice is the gauge field. Lattice gauge fields are SU(3) group elements in the fundamental representation, i.e. 3× 3 complex unitary matrices with determinant one. We associate one such element to each directed link on the lattice. If i and i+μ̂ label two neighboring lattice sites along the direction μ∈{1,2,3,4} the associated link is denoted by U_μ,i as in the figure.




    < g r a p h i c s >






 §.§ Lattice configurations


The collection of all lattice links is a lattice configuration. The basic linear operator in lattice QCD is the matrix of such a configuration along the direction μ:

    U_μ=[ U_μ,1            ; U_μ,2      ;           ⋱; U_μ,N ] ,

where N is the total number of lattice sites and the ordering is lexicographical. If N_1,N_2,N_3,N_4 are the number of sites along each direction, the total number of sites is N=N_1N_2N_3N_4. Note that U_μ is a block diagonal matrix with blocks of 3× 3 size.

QCDLAB follows the same data structure. For example, a random gauge field, which is appropriate for QCD at strong copuling can be generated by the following routine:



Note the sequential loop creating the gauge fields one by one. It is the only instance where QCDLAB uses such a loop in connection to degrees of freedom. This routine is called only once, usually at the begining of the simulation code. The user can avoid its call by simply starting with identity matrices and create random gauge fields using the simulation routine. This is controlled by the  function. Note also that the set of four gauge fields is stored in a single matrix .



 §.§ Shift operators


An important operator on the lattice is the permutation operator that shifts lattice sites along the positive direction μ:

    T_μ=[ 1; 1; ⋱; 1 ] .

In total we have four such operators, one for each direction. In four dimensions these are built using Kronecker products with identity matrices I_N_1,I_N_2,I_N_3,I_N_4:

    E_1   =   I_4⊗ I_3⊗ I_2⊗ T_1⊗ I_3
    
    E_2   =   I_4⊗ I_3⊗ T_2⊗ I_1⊗ I_3
    
    E_3   =   I_4⊗ T_3⊗ I_2⊗ I_1⊗ I_3
    
    E_4   =   T_4⊗ T_3⊗ I_2⊗ I_1⊗ I_3 .

Note that the extra Kronecker product with the identity 3× 3 matrix is neccessary in order to accomodate the space of gauge fields. The following routine creates the required operators:



Like in the case of gauge fields the set of four shift operators is stored in a single matrix .



 §.§ Wilson action


With the above operators we can write down the action of the SU(3) lattice theory as proposed by Wilson:

    S_gauge(U)=-1/g^2∑_μνtr (U_μ E_μ)(U_ν E_ν)(U_μ E_μ)^*(U_ν E_ν)^* ,

where the star symbolises the matrix Hermitian conjugation, the trace is taken in the 3N dimensional space and g is the bare coupling constant of the theory. In lattice gauge theory it is standart to use the inverse coupling constant:

    β=6/g^2 .

An important observable in QCD is the Wilson loop. The elementary Wilson loop, called plaquette, is the colorless product of gauge fields around an elementary suqare:

    P_μν,i=1/3Re tr_CU_μ,iU_ν,i+μ̂U_μ,i+ν̂^*U_ν,i^* ,

where the trace in performed in the color space.

(100,180)(0,0)
(150, 40)( 1, 0)100
(250, 40)( 0, 1)100
(250,140)(-1, 0)100
(150,140)( 0,-1)100

(150, 30)i
(250, 30)i+μ̂
(250,140)i+μ̂+ν̂
(120,140)i+ν̂

(190, 20)U_μ,i
(260, 80)U_ν,i+μ̂
(190,150)U_μ,i+ν̂^*
(120, 80)U_ν,i^*


It is straightforward to see that the product of matrices in the Wilson action can be written as a sum over all plaquettes: 

    S_gauge(U)=-β∑_i,μ>νP_i,μν .

The sum in the right hand side can be computed in QCDLAB by the following routine:



It follows directly the definition of the action in terms of sparse matrices, which makes the computation very efficient. Since the color trace is unnormalized in the routine the Wilson action is computed by calling .



 §.§ Dirac operator


Let D(U) be the Wilson formulation of the Dirac operator describing one quark flavor with bare mass m in the background gauge field configuration U:

    D(U)=(m+d)I-1/2∑_μ[(U_μ E_μ)⊗(1-γ_μ)+(U_μ E_μ)^*⊗(1+γ_μ)] ,

where d=4 and γ_μ are 4× 4 Dirac matrices obeying the Dirac-Clifford algebra in the Euclidean signature:

    {γ_μ,γ_ν}=δ_μν .

We assume also periodic boundary conditions in each direction. To see why it works we set U_μ=I and go to momentum space, in which case E_μ(p)=e^ip_μ and therefore:

    D(p)=m+i∑_μγ_μsin p_μ+∑_μ(1-cos p_μ) .

It is clear that for small momenta D(p)→ m+i+p^2/2+O(^3), whereas at other corners of the Brillouin zone there are 15 additional heavy flavors with masses m+2,m+4,m+6,m+8 and spin structure described by different sets of gamma-matrices. Therefore, at small momenta, Wilson fermions describe a single flavor of fermions and break chiral symmetry even at m=0 by the p^2/2 term. We will discuss this issue further in later sections. Here is the routine that implements Wilson fermions:






 §.§ A first algorithm


One special task in QCDLAB is the exponentiation of su(3) algebras. The concrete form of an su(3) algebra associated to a SU(3) gauge field in the fundamental representation is a 3× 3 anti-Hermitian traceless matrix. We have the following task: given a block diagonal matrix P_μ of order 3N with non zero su(3) algebra blocks we would like to compute the gauge field matrix:

    U_μ=e^P_μ

without using loops over the lattice sites. Here is an algorithm that completes this task:



It is an implemetation of the exponential power expansion:

    e^P_μ=∑_k=1^nP_μ^k/k!+O[P_μ^(n+1)]

truncated at order n=24 using the Horner algorithm. The order is chosen such that the resulting gauge fields are SU(3) matrices in the working precision of GNU Octave. If in doubt, the user should use the routine . There are more efficient implementations if we were to write the routine in C++. In this case one can exponentiate su(3) algebras one at a time using the algorithm behind the  function of the GNU Octave.



§ QCD PATH INTEGRAL


In this paper we focus in the simulation of pure Yang-Mills theory. Simulation of lattice QCD in this approximation, known as the quenched approximation, neglects screening coming from quark-antiquark pairs. It delievers very fast QCD properties such as the linear rising potential and hadron spectrum. Therefore, our task is the evaluation of the path integral:

    Z=∫∏_μ,idU_μ,i e^-S_gauge(U) ,

where dU_μ,i denotes the the SU(3) group integration measure, which is asummed to be gauge invariant. Its concrete form is unimportant in the algorithms used in QCDLAB.



 §.§ Hybrid Monte Carlo Algorithm


The HMC algorithm <cit.> starts by introducing su(3) conjugate momenta matrices P_μ to gauge fields. Gauge field configuations are generated by integrating classical field equations with Hamiltonian:

    H(P,U)=-1/4tr ∑_μP_μ^2-β/6∑_μνtr (U_μ E_μ)(U_ν E_ν)(U_μ E_μ)^*(U_ν E_ν)^* .

The extra one half in the normalization of the kinetic energy comes from the normalization of Gell-Mann matrices, which are adopted as su(3) algebra generators in the calculation of momentum matrices. The kinetic energy is computed by the following routine:



The first equation of motion is:

    U̇_μ=P_μU_μ .

Since H is an integral of motion, the second equation is derived by the equation:

    0=Ḣ=-1/4∑_μtr P_μṖ_μ
    -β/6∑_μν(U̇_μ E_μ)(U_ν E_ν)(U_μ E_μ)^*(U_ν E_ν)^*+h.c.   .

Substituting for U̇_μ the first equation of motion:

    0=Ḣ=-1/2∑_μtr P_μ{1/2Ṗ_μ+β/3∑_ν(≠μ)[ P_μν^(1)+ P_μν^(2)]}+h.c.

with P_μν^(1)=(U_μ E_μ)(U_ν E_ν)(U_μ E_μ)^*(U_ν E_ν)^* ,  P_μν^(2)=(U_μ E_μ)(U_ν E_ν)^*(U_μ E_μ)^*(U_ν E_ν) one gets the second equation of motion:

    1/2Ṗ_μ=-β/3∑_ν(≠μ)[ P_μν^(1)+ P_μν^(2)] .

Since P_μν matrices are 1× 1 loops around adjecent plaquettes that share a common link, they are block diagonal matrices of 3× 3 blocks. However, these blocks are not guaranted to be su(3) valued. Therefore, the force exerted at the gauge field U_μ,i is the traceless anti-Hermitian part of Ṗ_μ,i:

    F_μ,i=1/2(Ṗ_μ,i-Ṗ_μ,i^*)-1/3tr_C 1/2(Ṗ_μ,i-Ṗ_μ,i^*) .

The force is implemented in the following routine:




Having the equations of motion the next step is to build a trajectory using the leapfrog integration scheme:


    U_μ(t+Δ t/2)   =   e^P_μ(t)Δ t/2U_μ(t)
    
    P_μ(t+Δ t)   =   P_μ(t)+F_μ(t+Δ t/2)Δ t
    
    U_μ(t+Δ t)   =   e^P_μ(t+Δ t)Δ t/2U_μ(t+Δ t/2)


with a Δ t step size and a trajectory length τ. At t=0 su(3) momenta are taken to be Gaussian su(3) algebras: given eight independently distributed standard Gaussian variables at each lattice site and direction the routine  computes the corresponding momentum matrices. Note the half step updates of guage fields: it is expected that the force requires more flops than the exponentiation.

The algorithm ends by correcting for the non-conservation of the Hamiltonian using Metropolis et.al. with acceptance probability:

    P_acc({P(0),U(0)}→{P(τ),U(τ)})
    =min{1,e^-[H(τ)-H(0)]}

Upon rejection, one goes back to t=0 and refreshes momenta. This ends the description of the Hybrid Monte Carlo algorithm. One implicit and important assumption of QCDLAB is that the  function of GNU Octave suffices its purpose. The simulation routine of QCDLAB is:




Now we have enough programs to start exploring QCD. In the next section we begin with the string tension computation.



§ QCD STRING


In QCD, the energy between two static charges at large enough separation R follows the string law:

    V(R)=KR ,

where K is the string tension. On the lattice we measure the dimensionless string tension K̂=a^2K, where a is the lattice spacing. Using the string tension value from the Regge slopes, K=(440 MeV)^2, one can set the physical scale:

    a=197 MeV fm/440 MeV √(K̂) ,

where 197 MeV fm=1 is the energy-length conversion factor. Scale setting can be performed using any other physical quantity such as a hadron mass. QCD is compared to experiment by extrapolating dimensionless ratios of physical quantities at vanishing lattice spacing keeping the physical lattice size large enough to fit the physics. But how do we measure the string tension? We use an important lattice observable, the Wilson loop, which we deal with next.



 §.§ Wilson loop


We already know how to measure a 1× 1 Wilson loop, or the plaquette. Alghough one can measure all sorts of Wilson loops on the lattice, we restrict ourseleves to rectangle Wilson loops of dimensions R× T. Let V_1 and V_2 be the product of matrices along the R and T directions respectively:

    V_1=(U_1E_1)⋯(U_1E_1)_R times ,        V_2=(U_2E_2)⋯(U_1E_2)_T times .


(100,180)(0,0)
(150, 40)( 1, 0)100
(250, 40)( 0, 1)100
(250,140)(-1, 0)100
(150,140)( 0,-1)100

(120, 30)(0,0)
(250, 30)(R,0)
(250,140)(R,T)
(120,140)(0,T)

(190, 20)V_1
(260, 80)V_2
(190,150)V_1^*
(120, 80)V_2^*


Then, the (unnormalized) R× T Wilson loop is:

    W(R,T)=∑_μ≠νRe tr V_1V_2V_1^*V_2^* .

The foolowing routine is a direct implementation of the formula.






 §.§ Area law


The vacuum expectation value of R× T Wilson loop is the correlation function of the static quark-antiquark propagator separated with T lattice sites. Since, the large time behavior is dominated by the ground state contribution:

    W(R,T)=⟨ W(R,T)⟩≃ e^-V(R)T ,

the string behavior of the potential V(R)=KR is observed if the Wilson loop falls off exponentially with the area of the loop RT. A direct way to measure the string tension is by means of Creutz ratios:

    χ(R,T)=-logW(R+1,T+1)W(R,T)/W(R+1,T)W(R,T+1)

at large R and T values. Note that W's are sample averages and error propagation is not straightforward. A proper way to compute the error is using partial sample averages. This requires a large sample volume, which is often not available. An important trick, used as a short cut, is data resampling, or the so called bootstrap resampling. We will expain shortly, a bootstrap variant which is widely used in lattice QCD, the jackknife method.



 §.§ Jackknife resampling


Let us suppose we are given the real data vector x with volume n. The resampled jackknife data is the linear map:

    x^(J)=x̅-x-x̅/n-1 ,

which conserves the sample average x̅. If we require further the conservation of variance:

    Var[x^(J)]=Var(x)        ⇔        Cx^(J)-x̅^2=1/n-1x-x̅^2 ,

we should should choose the normalization factor C=n-1. As it is clear by inspection, the elements or x^(J) are partial sample averages of x:

    x_i^(J)=1/n-1∑_k(≠ i)x_k .

This way, we gain n-1 more sample averages of primary data than in the case without resampling. Thus, the physical quantity of interest, for example a Creutz ratio, is computed for all individual elements of x^(J) as opposed to a single sample average that was available originally. If y is the vector such estimations, its sample variance normalization is inherited from the corresponding variance of resampled data:

    Var(y)=(n-1)y-y̅^2 .

The QCDLAB routine that computes Creutz ratios is . Applying it to a small sample of 40 Wilson loops of maximal linear size 4, obtained on 6^3× 12 lattices, we get:

    χ(2,2)=0.372(3) ,        χ(3,3)=0.26(2) ,        χ(4,4)=0.2(1)

with the corresponding estimation of the lattice spacings, in fm units:

    a_(2,2)=0.273(1) ,        a_(3,3)=0.23(1) ,        a_(4,4)=0.19(7) .

Normally we should rely on the results of large Wilson loops. We see however, that the a_(4,4) value has a large error, so that the compromise is to select the a_(3,3) value as the estimation of our scale.



 §.§ Quark-antiquark potential


A standart way to measure the string tension on the lattice is measuring the quark-antiquark potential. In order to extract the potential from the Wilson loops one usually relies on effective potentials:

    V_eff(R,T)=-logW(R,T+1)/W(R,T) .

For fixed R we select as V(R) the median value of V_eff(R,T) over all T sizes of Wilson loops, whereas the corresponding error is computed using the jackknife method. These data are fitted to the general form:

    V(R)=V_o+α/R+KR ,

where V_o and α are two more constants in addition to the string tension K. Such a procedure is coded in the routine . Using the same set of Wilson loops as in the case of Creutz ratios the routine produces the 3-sigma band plot:




    < g r a p h i c s >




as well as the results:

    χ=0.23(2) ,            a_pot=0.22(1) fm .

From the plot we read that the potential is around 1.2 GeV if it is extrapolated at 1 fm separation. We observe also that the results overlap within 1-sigma with those obtained using Creutz ratios.



§ HADRON SPECTRUM


A basic computation in lattice QCD is the hadron spectrum. We will ilustrate it in the case of low lying mesons such as the pion and rho. Quark propagators are computed using the  routine. It calls the BiCGstab algorithm <cit.> as a Dirac solver:

    q=D^-1b ,

where D is the Wilson operator and b a point source at the origin of the lattice for each color and Dirac spin. Therefore, at each lattice site x the propagator q_x is a 12× 12 matrix. The pion and rho propagators are defined as:

    G_π,x=Tr q_xq_x^* ,        G_ρ,x=Tr γ_5γ_kq_xγ_kγ_5q_x^* ,

where the trace and Hermitian conjugation is performed in the tensor product space of color and spin. We sum over space-like lattice sites in order to get particle masses. For example, at long Euclidean time separation T the pion propagator is dominated by its ground state contribution:

    G_π,T=∑_x⃗G_π,(T,x⃗)≃ Ce^-m_π T .

Since we simulate with periodic boundary conditions in all directions the propagator decays exponentially also with respect to reflected times, which are translated by the lattice size N_4:

    G_π,T≃ C[e^-m_π T+e^-m_π(N_4-T)]∝cosh(T-N_4/2) .

Therefore, the routine  symmetrizes propagators with respect to the origin, which is actually at T=1:






 §.§ Effective masses


In complete analogy to the quark-antiquark potential we compute the effective masses of mesons as:

    M_eff(T)=-logG_π,T+1/G_π,T

and take the median value over all T values as the actual M_eff. The meson masses squared are then fitted against the bare quark masses using a linear model:

    M_eff^2=c_o+c_1m ,

where c_o,c_1 are unknown constants. With Wilson fermions we define the chiral limit at the vanishing pion mass. This procedure can be implemented by first calling the routine  with pion data to find the critical quark mass m_c. Then, the routine  is called with rho data and critical quark mass as input. It returns the lattice spacing using M_ρ=770 MeV. Finally, the routine  is called once more with pion data and the lattice spacing as input. Using the same configurations as before we get the 1-sigma band plot for pion and rho masses:




    < g r a p h i c s >




We have used a quadratic fit for the rho mass. The plot shows that the pion mass squared vanishes linearly with the quark mass in the chiral limit:

    M_π^2=c_1(m-m_c) .

The nonzero value of m_c is an artifact of Wilson fermions. For chirally symmetric fermions m_c should be zero. We these data, the estimated lattice spacing:

    a_ρ=0.23(1) fm

and the one estimated using the quark-antiquark potential overlap within 1-sigma errors.



§ AUTOCORRELATIONS


An important issue that one must take into the consideration in error reporting are autocorrelations in the Monte Carlo time series. If x^(1)=( O_1, O_2,…, O_n) is the time series vector of an observable O one measures the autocorrelation function between x^(1) and the time forwarded samples x^(1),x^(2),…,x^(t):

    f_j, O=C∑_k=1^n-t+1(x_k^(1)-x̅^(1))(x_k^(j)-x̅^(j)) ,        j=1,2,…,t

as well as the integrated autocorrelation time τ_int, O <cit.>:

    τ_int, O=1/2+∑_j=2^n(1-j-1/n)f_j, O .

The normalization constant C is chosen such that f_1=1. The right hand side may be approximated by the sum:

    τ_int, O≈1/2+∑_j=2^tf_j, O

assuming that the data volume is much larger than the cutoff t. QCDLAB routine that computes autocorrelations is :



A proper estimation of autocorrelations should also ensure that the integrated autocorrelation times are small compared to t, and the corresponding error is computed on a large data set. In our simulation example, plaquette decorrelates in 5(2) Hybrid Monte Carlo trajectories and saved configurations are separated by 100 trajectories.



§ MINIMALLY DOUBLED CHIRAL FERMIONS


In this section we describe new features included in the version 2.1 of QCDLAB. When studying the spontaneous breaking of chiral symmetry with Wilson fermions we encounter a nonzero critical bare quark mass, which is a consequence of the explicit breaking of the chiral symmetry of the Wilson-Dirac operator. One way to circumvent this problem si by using Ginsparg-Wilson fermions which posses exact chiral symmetry on the lattice <cit.>. Such fermions, which are part of the QCDLAB versions 1.0 and 1.1  <cit.>, are compuationally complex and we skip them in the follwing. Instead, we use here chiral fermions with broken hypercubic symmetry <cit.>. These fermions describe a degenerate isospin doublet and are ideal for studying QCD with u and d quarks. However, due to broken hypercubic symmetry they lead to counterterms in the action <cit.>. We use the version of reference <cit.> since it has a negligibly small gauge action counterterm, which we set to zero as a first approximation. The structure of the Dirac operator resembles the one of the Wilson operator:

    D_BC(U)=mI+i(-2+c_3)Γ+1/2∑_μ[(U_μ E_μ)⊗(γ_μ+iγ_μ')+(U_μ E_μ)^*⊗(-γ_μ+iγ_μ')] ,

where the new set of gamma matrices γ_μ'=Γγ_μΓ is associated to the doublet partner fermion and Γ=∑_μγ_μ/2=∑_μγ_μ'/2. Apart from a negligibly small counterterm, which we have set to zero, the counterterm ic_3Γ is used to restore the hypercubic symmetry of the theory. In the follwoing figure we see that the neutral pion mass, the only Goldstone boson of the theory, flattens at c_3=0.2 (we have fixed the bare quark mass at the small value 0.01):




    < g r a p h i c s >




We have used the same meson operators as with Wilson fermions keeping in mind the two-fold degeneracy of the isospin doublet fermion. For example, if the fermion doublet ψ is written in terms of flavor singlet fields u and Γ d, i.e. ψ=u+Γ d, the pion operator admits the expression of the neutral pion operator:

    ψ̅γ_5ψ=u̅γ_5 u+d̅ Γγ_5Γ d=u̅γ_5 u-d̅γ_5d

(Γ^2=1 and it anticommutes with γ_5). Note that cross terms vanish identically since u and d are defined on the orthogonal subspaces of the flavor doublet field ψ. With the value of c_3 being fixed at 0.2 one may compute the pion and rho masses for a range of bare quark masses. Below we plot the fitted data with the 1-sigma band:




    < g r a p h i c s >




Note that we have included a quenched chiral logarithm in the pion mass model:

    M_π^2=c_1 m+c_2 mln m ,

since it gives a better fit of the data. The lattice spacing value a_ρ=0.20(3) fm overlaps with the value found using Wilson fermions within the 1-sigma error.

In summary, we have described briefly the main routines of QCDLAB, versions 2.0 and 2.1, as well as its use.



 §.§ Acknowledgement


The author thanks Travis Whyte who found a missmatch of shift operators at the  directory and inside the  function.

10

gross_wilczek D. J. Gross, Frank Wilczek, Ultraviolet Behavior of Nonabelian Gauge Theories, Phys. Rev. Lett. 30 (1973) 1343.

politzer H. D. Politzer, Reliable Perturbative Results for Strong Interactions?, Phys. Rev. Lett. 30 (1973) 1346. 

wilson K. G. Wilson, Confinement of Quarks, Phys. Rev. D10 (1974) 2445.

kogut_susskind J. B. Kogut, L. Susskind, Hamiltonian Formulation of Wilson's Lattice Gauge Theories, Phys. Rev. D11 (1975) 395.

creutz M. Creutz, Confinement and the Critical Dimensionality of Space-Time, Phys. Rev. Lett. 43 (1979) 553-556.

octave .

hmc S. Duane et. al., Hybrid Monte Carlo, Phys. Lett. B195 (1987) 216.

van_der_vorst H. A. Van der Vorst, Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems, SIAM J. Sci. and Stat. Comput. 13(2) (1992), 631.

gutknecht M. H. Gutknecht, Variants of BICGSTAB for Matrices with Complex Spectrum, SIAM J. Sci. Comput. 14(5) (1993) 1020.

sokal A . D. Sokal, Bosonic Algorthms, in Quantum Fields on the Computer, M. Creutz, editor, World Scientific (1992) 211.

neuberger H. Neuberger, Exactly massless quarks on the lattice, Phys.Lett. B417 (1998) 141.

qcdlab1 A. Boriçi, QCDLAB: Designing Lattice QCD Algorithms with MATLAB, , A. Boriçi, Speeding up Domain Wall Fermion Algorithms using QCDLAB, .

kw L.H. Karsten, Lattice Fermions in Euclidean Space-time , Phys. Lett. B 104 (1981) 315; F. Wilczek, On Lattice Fermions, Phys. Rev. Lett. 59 (1987) 2397.

bc M. Creutz, Four-dimensional graphene and chiral fermions , JHEP 0804 (2008) 017; A. Boriçi, Creutz fermions on an orthogonal lattice, Phys. Rev. D. 78 (2008) 074504.

capitani_et_al S. Capitani, M. Creutz, J. Weber, H. Wittig, Renormalization of minimally doubled fermions, JHEP 1009 (2010) 027.



