# Sum on a lattice where variables can only take two values (Ising lattice gauge theory)

I want to study the Ising Gauge Theory which has the Hamiltonian $$$$\mathcal{H} = -J \sum_{n,\mu,\nu} \sigma_3(n,\mu)\sigma_3(n+\mu,\nu)\sigma_3(n+\mu+\nu, -\mu)\sigma_3(n+\nu, -\nu)$$$$ where $$n$$ labels a node and $$\mu, \nu$$ are unit vectors along the lattice. The spins ($$\sigma_3 = \pm 1$$) are placed on the bonds between two nodes. I wanted to study this model in 2 and 3 dimensions by calculating it's partition function where I naively want to define the $$i$$th plaquette as

plaq[i_] = sxb[i] syb[i] sxt[i] syt[i]


where sxt[i] is the spin on the x axis for the 'top' bond of the $$i$$th square and sxb[i] is the spin on the bottom bond. And similar definition for the y bonds. But this term cannot be summed easily over the squares since Sum does not sum over both the spin variables and their label i. Hard coding the sum is a very inefficient way to do this, so I was looking for better ways. Is there a way Mathematica can do this sum? Efficiency is not an issue since I will not examine very large lattices. Also, help for the same with a 3D lattice would be greatly appreciated. Thanks a lot.

One option is to specify each node as a pair of integers, {x,y}. You might also consider placing your spins at the nodes rather than on the edges. There will be four spins/plaquette either way, and the spin associated with a link can be though of as the product of the two spins on the nodes.

Start by randomly assigning spins to your lattice. I'll assume it runs from {1,xMax} and {1,yMax}:

lattice = 2 RandomInteger[{0,1},{xMax,yMax}] -1/2;


The generation is handled this way to omit the zero that results if I simply try to use RandomInteger[{-1,1}].

You can define a plaquette by passing the coordinates of the lower left corner, with spins positioned at the four corners. To retrieve the current values, define getPlaquette:

getPlaquette[{x_,y_}] := lattice[[Sequence@@#]]& /@{{x,y},{x+1,y},{x+1,y+1},{x,y+1}};


To set an individual node's spin to a new value, use setNodeSpin:

setNodeSpin[{x_,y_}, value_Integer] := lattice[[x,y]] = value;


You can also write a function to set all four spins at the corner of a plaquette using by calling setNodeSpin multiple times with the appropriate coordinates passed in.

Your Hamiltonian sum then becomes a sum over {x,1,xMax}, {y,1,yMax}.

If you need to specify a single identifier for each node, you can take the {x,y} pair and convert it in a similar manner to that used to index a multi-dimensional array as a 1D array, for example:

nodeId[x_,y_] := (x-1)*xMax + y;

• Getting {1, 1, 1, -1} on the bonds cannot be done by multiplying spins on the nodes. Nov 18, 2023 at 12:35