How to write a naive Dirac matrix for Lattice QCD?

Question

I'm trying to write down the naive Dirac matrix (with fermion doubling) for a LQCD simulation with one quark, for now. I initialized the $SU(3)$ gauge field and the quark field. The quak field has 4 space-time indices, 1 index for spinor components (from 1 to 4) and 1 index for the color (from 1 to 3), if you agree. My question is: how to couple together all these indices in Dirac matrix? And how many indices does the Dirac matrix have? I'm thinking of something like this (in Python):

spinor=4 
mu=4 #(all 4 directions in space-time) 
color=3 
su3=3 

quark=np.zeros((Nx, Ny, Nz, Nt, spinor, color), complex) 
gaugeSU3=np.zeros((Nx, Ny, Nz, Nt, mu, su3, su3), complex)

Once filled all quark and gauge fields with appropriate values (for the gauge there are 4 SU(3) matrices for each space-time point), how to treat all the indices in Dirac matrix? Simulations of SU(3) pure gauge theory I performed are coherent with the literature, so I suppose it works, but how to couple the field with quarks? I add here the code I wrote for naive fermions; It seems to be incorrect because the Dirac matrix isn't a sparse matrix but all elements are filled with values different from zero:

 def DiracMatrix(U, psi, D):
   #here U is the gaugeSU3, psi is the quark field, D the Dirac 
   # matrix, Dirac is the spinor indices
   m = 0.2

    for x in range(Nx):
      for y in range(Ny):
         for z in range(Nz):
            for t in range(Nt):
                for alpha in range(Dirac):
                   for beta in range(Dirac):
                       for a in range(color):
                           for b in range(color):
                               for mu in range(4):
                                    a_mu = [0, 0, 0, 0]
                                    a_mu[mu] = 1
                                    D[x, y, z, t, alpha, beta, a, b] += 0.5 * (
                                        gamma[mu][alpha, beta]
                                        * U[x, y, z, t, mu, a, b]
                                        * psi[
                                            (x + a_mu[0]) % Nx,
                                            (y + a_mu[1]) % Ny,
                                            (z + a_mu[2]) % Nz,
                                            (t + a_mu[3]) % Nt,
                                            alpha,
                                            a,
                                        ]
                                        - U[
                                            (x - a_mu[0]) % Nx,
                                            (y - a_mu[1]) % Ny,
                                            (z - a_mu[2]) % Nz,
                                            (t - a_mu[3]) % Nt,
                                            mu,
                                            a,
                                            b,
                                        ]
                                        .conj()
                                        .T
                                        * psi[
                                            (x - a_mu[0]) % Nx,
                                            (y - a_mu[1]) % Ny,
                                            (z - a_mu[2]) % Nz,
                                            (t - a_mu[3]) % Nt,
                                            beta,
                                            b,
                                        ]
                                        + m
                                    )

return D

Please help me to comprise better this step! Thanks

Answer 1

On a 4d cubic lattice, the naive Dirac action just has the usual 4-by-4 $\gamma_0$ matrix on the time-direction links and the 4-by-4 $\gamma_1$ on the $x$-direction links, etc. The gamma's talk only to the spinor indices.

At each site you will store 4 spin compoments. If examine your code you will see that there are four sets of variables, variables in the same set talk to each other through the links, but varibles in different sets do not talk to each other, so you need keep only one per site. If you do this you end up with Kogut-Sussing "staggered fermions"

You need separate ${\rm SU}(3)$ colour matrices $U_\mu$ for the gauge coupling on the links in the $\mu$ direction. Do you not have a book that explains all this?