I want to study the Ising Gauge Theory which has the Hamiltonian \begin{equation} \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) \end{equation} 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.


1 Answer 1


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;
  • $\begingroup$ Getting {1, 1, 1, -1} on the bonds cannot be done by multiplying spins on the nodes. $\endgroup$
    – QFTheorist
    Nov 18, 2023 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.