モノポール？

元の論文は説明を引用文献に任せてあるんで単体で読んでも分かりにくい。こちらを読む必要ありそう
G.E. Volovik, J. Phys. C: Solid State Phys. 20, L83-L87 (1987).

上記論文の自分の解釈ですが、Eq.1は自由電子のAction＋ローカルなスピン偏極の４乗。Φはその場所そのスピンの生成消滅演算子。
αはたとえば=↑、↓の二つ。最後の項がなければ相互作用なくてそれぞれ平面波が解。最後の項のためにスピン偏極≠０になる。空間的にその向きが変わるけど大きさは同じとしてみる。それを真空として励起を考える。
最後の項は定数になって消えるけど場所ごとに基底が変わるから前の２つも変わってEq.(4)になる。 前２つだけ見ればこれはまさにベクトルポテンシャルA中の自由電子。 でもこのAは本当の電磁場じゃないのでスピン偏極の時間空間配置をいじればdivB≠０にできる。

JPSJ論文の方は場所ごとに違う方向に磁場がかかってて、そっちむきのスピンのエネルギーを下げる。座標を場所ごとに取り直して磁場の方向を基準にする。
普通は波動関数とか隣と同じ値になろうとする。でも上みたいなことをするとひねりが加わって、隣と若干違った値になろうとする。具体的にはスピン↑と↓の二成分あるのを、SU(2)行列でひねったものに揃おうとする。そういう隣同士にひねりを加えるようなのをゲージ場という。
具体的にはexp(iパウリ行列)になる。∂μを∂μ-i{{Az,Ax+iAz},{Ax-iAy,-Az}}で置き換える。Ax,Ay,Azは方向μごとに違うので添え字μもつく。

とりあえずストレートに時間発展を計算するコード書いてみた。HはA+B、それぞれ交換しない、と表せるけど、単純にexp(idt A)exp(idt B)をかけていってる。基底はゲージ変換とか考えずにもとのまま。
はじめからHを厳密に対角化すれば一番効率的ですが、まあそれはそれ。

// Solve d |phi1,phi2>/dt = i H |phi1,phi2>
// H= Laplacian + local magnetic field described by Pauli matrixes

// define this to use imaginary time update for ground state calclation
//#define DAMP

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//system size
#define L 32
// 2D system
#define VOL (L*L)

typedef double complex[2];
// wavefunction for spin up/down
complex phi1[VOL],phi2[VOL];
// 4 neighbor index
int nnidx[VOL][4];
// local field
double hfield[VOL][3];

typedef complex c2matrix[2][2];
// local spin rotaion matrix : exp(i dt h_a sigma_a)
c2matrix umat[VOL], pmat[VOL];

#define XY2I(x,y) ((y)*L+(x))
#define czero(c) {c[0]=c[1]=0;}
#define cset(c,R,I) {c[0]=R;c[1]=I;}
// c=ab
void cmul(complex a, complex b, complex c)
{
   c[0]=a[0]*b[0]-a[1]*b[1];
   c[1]=a[0]*b[1]+a[1]*b[0];
}
//a+=b
void cinc(complex a, complex b)
{
    a[0]+=b[0];
    a[1]+=b[1];
}
// cij = sum_k aik bkj
void mmul(c2matrix a, c2matrix b, c2matrix c)
{
   int i,j,k,l;
   complex cc;
   for(i=0;i<2;i++){
   for(j=0;j<2;j++){
       czero(c[i][j]);
       for(k=0;k<2;k++)
       {
           cmul(a[i][k], b[k][j], cc);
           cinc(c[i][j], cc);
       }
   }}//ij
}
// u2=a.u
void mvmul(c2matrix a, complex u[2], complex u2[2])
{
   int i,j,k,l;
   complex cc;
   for(i=0;i<2;i++){
       czero(u2[i]);
       for(k=0;k<2;k++)
       {
           cmul(a[i][k], u[k], cc);
           cinc(u2[i],cc);
       }
   }//ij
}

// spin rotation matrix  um=exp(i dt m_a sigma_a)
void calU(double mx, double my, double mz, double dt,  c2matrix um, c2matrix pm)
{
   double m0 = sqrt(mx*mx+my*my+mz*mz);
   double s=mz/m0;
   double cr=mx/m0;
   double ci=my/m0;
   double c=sqrt(mx*mx+my*my)/m0;


   c2matrix u,tu, tmp1,tmp2,dd;
   double r1=sqrt(c*c+(1-s)*(1-s));
   double r2=sqrt(c*c+(-1-s)*(-1-s));

   // spin Hamiltonian
   cset(pm[0][0],  mz*dt, 0);
   cset(pm[1][1], -mz*dt, 0);
   cset(pm[0][1], mx*dt, my*dt);
   cset(pm[1][0], mx*dt, -my*dt); 

   // c , (1-s)/e | c e,-1-s
   if(r1==0.0)
   {
       cset(um[0][0], cos(dt*mz),  sin(dt*mz));
       cset(um[1][1], cos(dt*mz), -sin(dt*mz));
       czero(um[0][1]);
       czero(um[1][0]);
       return;
   }
   cset(u[0][0], mx/r1, my/r1);
   cset(u[0][1], mx/r2, my/r2);
   cset(u[1][0], (1-mz)/r1, 0);
   cset(u[1][1], (-1-mz)/r2,0);

   int i,j,k;
   // transpose, bar
   for(i=0;i<2;i++){
   for(j=0;j<2;j++){
       tu[i][j][0]=u[j][i][0];
       tu[i][j][1]=-u[j][i][1];
   }}

   // diagonalized matrix
   for(i=0;i<2;i++){
   for(j=0;j<2;j++){
       czero(dd[i][j]);
   }}
   cset(dd[0][0],cos(dt), sin(dt));
   cset(dd[1][1],cos(dt),-sin(dt));
   // um=u.dd.tu
   mmul(dd,tu, tmp1);
   mmul(u,tmp1, um);
}

// init neighbor index
void innidx()
{
   int x,y;
   for(x=0;x<L;x++){
   for(y=0;y<L;y++){
       int x1=(x-1+L)%L;
       int x2=(x+1+L)%L;
       int y1=(y-1+L)%L;
       int y2=(y+1+L)%L;
       int i0=XY2I(x,y);
       nnidx[i0][0]=XY2I(x1,y);
       nnidx[i0][1]=XY2I(x2,y);
       nnidx[i0][2]=XY2I(x,y1);
       nnidx[i0][3]=XY2I(x,y2);
   }}
}
// Hphi = i dt H phi, H= Laplacian
void H_phi(complex phi[VOL], complex Hphi[VOL], double dt)
{
   int i,j,k;
   for(i=0;i<VOL;i++)
   {
       complex c;
       for(j=0;j<2;j++)
       {
           c[j] = -4*phi[i][j];
           for(k=0;k<4;k++) c[j] += phi[nnidx[i][k]][j];
       }
#ifdef DAMP
       Hphi[i][0] =  dt* c[0];
       Hphi[i][1] =  dt* c[1];
#else
       // mult i
       Hphi[i][0] =  dt* c[1];
       Hphi[i][1] = -dt* c[0];
#endif
   }
}
complex dphi[VOL];
complex phin[VOL];
complex phinn[VOL];

// phi += \sum_n=1^nmax (i dt H)^n/n!
void upd1(complex phi[VOL], double dt)
{
   int i,j,k,n;
   int nmax=10;

   for(i=0;i<VOL;i++) dphi[i][0]=dphi[i][1]=0;

   H_phi(phi, phin, dt);
#ifdef DAMP
   for(i=0;i<VOL;i++)
           for(j=0;j<2;j++) phi[i][j] += phin[i][j];
#else
   for(n=1;n<=nmax;n++)
   {
       for(i=0;i<VOL;i++){
       for(j=0;j<2;j++){
           phin[i][j]/=n;
           dphi[i][j] += phin[i][j]; // +(i dt)^n H^n|phi>/n!
       }}

       // phin = i dt H phin
       H_phi(phin, phinn, dt);
       for(i=0;i<VOL;i++)
           for(j=0;j<2;j++) phin[i][j] = phinn[i][j];
   }

   for(i=0;i<VOL;i++)
       for(j=0;j<2;j++) phi[i][j] += dphi[i][j];

#endif
}

// output for gnuplot (3d splot)
void dump()
{
   int i,j,x,y;
   double psum=0;
   for(x=0;x<L;x++){
   for(y=0;y<L;y++){
       int ix=XY2I(x,y);
       double c1r=phi1[ix][0];
       double c1i=phi1[ix][1];
       double c2r=phi2[ix][0];
       double c2i=phi2[ix][1];
       psum += c1r*c1r+c1i*c1i;
       psum += c2r*c2r+c2i*c2i;
       double r1=sqrt(c1r*c1r+c1i*c1i);
       double r2=sqrt(c2r*c2r+c2i*c2i);

       double mz= r1*r1-r2*r2;
       double mx= 2*(c1r*c2r-c1i*c2i);
       double my= 2*(c1r*c2i+c2r*c1i);
       double rr=r1*r1+r2*r2;

       double *h = hfield[ix];
       double th1=atan2(c1r,c1i);
       double th2=atan2(c2r,c2i);
       printf("%d %d %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
               x,y, rr,th1,th2,
		c1r,c1i,c2r,c2i,
	       mx,my,mz,
               h[0],h[1],h[2]);
   }
   printf("\n");
   }
   printf("#sum %f\n",psum/VOL);
}

void init_field()
{
   int x,y;
   for(x=0;x<L;x++){
   for(y=0;y<L;y++){
       int ix=XY2I(x,y);
       double dx=x-L/2;
       double dy=y-L/2;
       double dd=fabs(dx);
       if(fabs(dy)>dd)dd=fabs(dy);

       if(x==L/2 && y==L/2)
       {
           hfield[ix][2]=1;
           hfield[ix][0]=0;
           hfield[ix][0]=0;
           continue;
       }
       double th = atan2(dy,dx)*1;
       double wid=L/4;
       double ph = M_PI*(1-exp(-(dx*dx+dy*dy)*1.0/wid/wid));
       //double ph=dd/(L/2)*M_PI;
       hfield[ix][2]=cos(ph);
       hfield[ix][0]=sin(ph)*cos(th);
       hfield[ix][1]=sin(ph)*sin(th);
   }}
}

void init_u(double dt)
{
   int i;
   for(i=0;i<VOL;i++)
   {
       double *h = hfield[i];
       calU(h[0],h[1],h[2], dt, umat[i], pmat[i]);
   }
}

//rotate spin  |phi>'=exp(i dt h_a sigma_a) 
void upd2()
{
   int i;
   for(i=0;i<VOL;i++)
   {
       double p1[2][2];
       double p2[2][2];
       p1[0][0]=phi1[i][0];
       p1[0][1]=phi1[i][1];
       p1[1][0]=phi2[i][0];
       p1[1][1]=phi2[i][1];
#ifdef DAMP
       mvmul(pmat[i], p1, p2);
       phi1[i][0] -=p2[0][0];
       phi1[i][1] -=p2[0][1];
       phi2[i][0] -=p2[1][0];
       phi2[i][1] -=p2[1][1];
#else
       mvmul(umat[i], p1, p2);
       phi1[i][0]=p2[0][0];
       phi1[i][1]=p2[0][1];
       phi2[i][0]=p2[1][0];
       phi2[i][1]=p2[1][1];
#endif

   }
}
void normalize()
{
   double psum=0;
   int i;
   for(i=0;i<VOL;i++)
   {
      double *d1=phi1[i];
      psum+= d1[0]*d1[0] + d1[1]*d1[1];
      d1=phi2[i];
      psum+= d1[0]*d1[0] + d1[1]*d1[1];
   }
   psum=sqrt(VOL/psum);
   for(i=0;i<VOL;i++)
   {
      double *d1=phi1[i];
      d1[0]*=psum; d1[1]*=psum;
      d1=phi2[i];
      d1[0]*=psum; d1[1]*=psum;
   }
}

int main(int argc, char **argv)
{
   double h=atof(argv[1]);
   double dt=1e-2;
   int i;
   innidx();
   for(i=0;i<VOL;i++) {phi1[i][0]=1.0; phi1[i][1]=0;}
   for(i=0;i<VOL;i++) phi2[i][0]=phi2[i][1]=0;
 init_field();
 init_u(dt*h);

   for(i=0;i<5000;i++)
   {
       //Suzuki-Trotter 2nd order exp(dt B/2)exp(dt A)exp(dt B/2)
       // is equivalent to exp(dt B) exp(dt A) exp(dt B) exp(dt A) ...

       // Plane wave
       upd1(phi1,dt);
       upd1(phi2,dt);
       //spin
       upd2();
#ifdef DAMP
       normalize();
#endif
   }

   dump();
}