Maxwell方程式の表現論

Maxwell方程式の表現論(0)
http://d.hatena.ne.jp/m-a-o/20150125#p2
からの続き(本編)。

2021年追記)Cartan部分環の取り方が正しくないので、以下の結果は、ダメっぽい。正しいスカラー場の計算は、以下を参照
so(n+1,2)の極小表現の実現
https://vertexoperator.github.io/2021/05/05/hermite_BDI_minrep.html


Poincare代数の離散ヘリシティを持つmassless既約ユニタリ表現は、共形代数$su(2,2)$のladder representationに一意に拡張できて、twistor理論では、Dolbeaultコホモロジーを使って、この表現空間の実現が与えられた。一方で、共形変換はMaxwell方程式の解空間にも自然に作用していて、DolbeaultコホモロジーからMaxwell方程式の(適当なクラスの)解空間(正確には、自己双対Maxwell方程式と反自己双対Maxwell方程式それぞれの解空間)への表現の同値を与えるintertwinning operatorは、Penrose変換によって具体的に実現される。と言いたいのだけど、Penrose変換の像がどうなっているか、明示的に書いてある文献を見つけることができなかった。


というわけで、真空中の自己双対Maxwell方程式(SD-Maxwell方程式)と反自己双対Maxwell方程式(ASD-Maxwell方程式)の解の中で、ladder representationと同値な表現を与えるような基底を見つけることを考える。twistor空間上の"関数"は、物理的な解釈がしづらいものなので、Minkowski空間上の場に対して、ladder representationが実現できてしまえば、Penrose変換とかコホモロジーとか難しいものは忘れてしまってもよい



方針。$su(2,2)$のladder representationは、lowest weight representationであって、この既約ユニタリー表現を、純表現論的に記述することは(原理的には)難しくない。Verma加群の商としての直接的な記述は、例えば
THE MASSLESS REPRESENTATIONS OF THE CONFORMAL QUANTUM ALGEBRA
http://iopscience.iop.org/0305-4470/27/14/012
などで見ることができる。表現の構造は、これで分かったとして、SD/ASD-Maxwell方程式の解の中で、lowest weight vectorに対応するものを発見して、それに対して、無限小共形変換を、どんどん作用させていけば、ASD/SD-Maxwell方程式の解空間の基底が得られる(はず)。これがladder representationと同値であることを確認するには、lowest weight vectorが満たす最低ウェイト条件の確認と、singular vectorがちゃんと消えることを言えばいい。


lowest weight vectorは適当なPenrose変換の像の中にあると信じて探すのだが、
On the density of twistor elementary states.
http://projecteuclid.org/euclid.pjm/1102637077
に、elementary statesとは$K$-finite vector($K$は$SU(2,2)$の極大コンパクト部分群)と同じものであると書いてある。なので、elementary statesの像の中で探せばよい。勿論、最低ウェイト条件もsingular vectorが消えるという条件も線形微分方程式になるので、直接解いてもいいけど


#これらの基底は、いずれも"ハミルトニアン"(ここでは、Hamiltonianは、Poincare代数の時間並進演算子のこととする)の固有状態ではない(多分、エネルギー固有状態は、一つも存在しない。これは、量子力学に於いて、平面波は、ラプラシアンの"固有状態"であるが、二乗可積分でない、というのと似た事態だと思う)。ところで、Weinbergの教科書「場の量子論」1巻の2.5節の冒頭を見ると、"4元エネルギー・運動量テンソルは互いに可換だから、物理的状態ベクトルをこの4元ベクトルの固有値を使って表すことは自然だ."と書いてある。これは、連続スペクトルと固有値を同じように扱おうということだと思うけど、Weinbergの教科書に限らず、物理の教科書は、大体こういう議論で始まる気がする(Wignerの論文自体で、こういう議論がされていたっぽい。cf: arXiv:0809.4942)。


準備。まず、spinor場への共形変換の作用を具体的に知る必要がある。これは
Finite-component field representations of the conformal group
http://www.sciencedirect.com/science/article/pii/0003491669902784
で調べられている。論文では、spinor場とは限らない、一般の多成分場を扱っていて、共形変換の多成分場への表現は、適当な1-cocycleでtwistした誘導表現で実現されると考え、このcocycleが並進に対しては自明である時、スケール変換・Lorentz変換・特殊共形変換のなす群(の被覆群)の有限次元表現とcocycleが対応する。という具合であるが、今考えるのはspinor場であるので、特殊共形変換の作用は自明であるとして、結局、$Spin(3,1)$の有限次元表現から、cocycleを作り、誘導表現を得ることになる。結局、これは多くの場の理論の教科書に書いてあることと同じで、massless spinor場の場合、$Spin(3,1)$の有限次元表現として、$D(j,0)$あるいは$D(0,j)$の形のものを取る($j$は非負の整数or半整数で、それぞれhelicityが+jあるいは-jの場に対応する)。こうして、spinor場への共形変換の作用が分かる。$Spin(3,1)$の有限次元表現と$so(3,1)$の有限次元表現は、1:1に対応するので、以下、両者を混同する


#Poincare代数の作用は、Klein-Gordon作用素と可換(Klein-Gordon作用素はPoincare代数のCasimir元)であるけども、特殊共形変換やスケール変換は、Klein-Gordon作用素と可換ではない(けど、Klein-Gordon方程式の解を別の解を移す)ので、同じ"対称性"といっても、Poincare代数の場合と、共形代数の場合では、少し意味が違う(対称性という言葉は、濫用されすぎな気もする)


練習。"spin 1"(helicity +1/-1)の場合を考える前に、masslessなスカラー粒子の場合を眺める。この場合、spinor場の満たす方程式は、質量0のKlein-Gordon方程式。Penrose変換を計算する必要があるのだけど、幸いなことに
(Pre-)Hilbert spaces in twistor quantization
http://arxiv.org/abs/1210.0349
で必要な計算をしてくれているので、カンニングすることにする。結論を先に書くと、
\phi = \frac{1}{x_0^2 - x_1^2 - x_2^2 - x_3^2}
が最低ウェイトベクトルを与える。


まず、最低ウェイト条件+特異ベクトルが消える条件をChevalley生成元を使って書く。共形代数は、複素化すると、sl(4)ないしso(6)と同型なLie代数となり、これはrankが3である。最低ウェイト条件は、まず
f_1 \phi = f_2 \phi = f_3 \phi =0
で、Cartan部分環の作用は、非負の半整数(or整数)の組$(j_1,j_2)$と共形次元$d = j_1+j_2+1$を使って
(h_1 + 2 j_1)\phi = (h_2 -d - j_1 - j_2)\phi = (h_3 + 2 j_2)\phi= 0
と書ける(ことが分かっている)。スカラー粒子の場合は、j1=j2=0で、
h_1 \phi = (h_2 - 1)\phi = h_3 \phi = 0
となる。また、2つの特異ベクトルが存在し、それは、$e_1$と$e_3$をそれぞれ一回ずつ作用させて得られる。これが0でないといけないので
e_1 \phi = e_3 \phi = 0
という2つの条件も確認する必要がある。


#[宿題]ユニタリー性を満たすためには、j1=0かj2=0でないといけないが、ユニタリー性を諦めれば、j1とj2が共に0でないことも可能である。そして、そのような表現は、$so(3,1)$の有限次元既約表現$D(j1,j2)$に付随するspinor場と関連すると予想するのは自然に思える。例えば、ベクトルポテンシャル(あるいは、一般のベクトル場)は、表現$D(1/2,1/2)$に付随するspinor場である(とされている)。この場合、Gupta-Bleuler形式を思い出して、最低ウェイトVerma表現から、物理でよくやるように、ゴーストを取り除いて商表現を作ると、2つのladder表現の直和になるだろうか(特異ベクトルを除くのも、ゴーストを除くのも、ユニタリーになるように商を取るだけなので、作業としては同じようなもんである)。このVerma表現を純表現論的に記述することは、多分そんなに難しくない。一方で、電磁ポテンシャルの中で、(j1,j2)=(1/2,1/2)の最低ウェイトベクトル条件を満たすものを見出すことも、多分できる。この"最低ウェイトベクトル場"に共形代数を作用させて得られる表現は、Verma表現と同値なのか、それとも、その何らかの商表現となっているのだろうか。


次に、Chevalley生成元と、共形代数の"幾何学的な生成元"の関係を知らないといけない。天下り的であるが
e_1 = \frac{1}{2}(M_1 + i M_2 + i L_1 + L_2)
e_2 = \frac{1}{2}(-P_{0} - P_{3})
e_3 = -\frac{1}{2}(M_1 - i M_2 -i L_1 - L_2)
h_1 = L_3 - i M_3
h_2 = -(D + L_3)
h_3 = L_3 + i M_3
f_1 = \frac{1}{2}(-M_1 + i M_2 - i L_1 - L_2)
f_2 = \frac{1}{2}(K_{0} - K_{3})
f_3 = \frac{1}{2}(M_1 + i M_2 - i L_1 + L_2)
が答え。$D$,$P_a$,$K_a$,$M_i$,$L_i$はそれぞれ、スケール変換、並進、特殊共形変換、空間回転、Lorentz boostのgenerator。そして、これらの生成元のspinor場への作用は、既に分かっているので、最低ウェイト条件+特異ベクトルの消滅を確認するのは、ただの計算練習となる。


(注意)スケール変換と特殊共形変換の生成子には、cocycleを作る時に使ったスケール変換の一次元表現から来る"余分な"項δ(論文"Finite-component field representations of the conformal group"参照)が付いている。δ=0とすると、普通の特殊共形変換の無限小生成子が得られるけど、スカラー粒子の場合でも、δは0でない(0としてしまうと、lowest weight vectorは特殊共形変換で不変でなくなり、困ったことになってしまう。当初スカラーの場合はδ=0でいいと思いこんでいて、半日くらい悩んでしまった。。)


本題。スカラーの場合が分かれば、"spin 1"の場合も、殆ど同じようなものである。spin 1の場合は、正ヘリシティ・負ヘリシティのそれぞれの場は、3成分を持つので、計算は面倒になる。面倒になるけど、計算自体は高校生でもできる程度のものである。末尾に、Risa/Asirによる計算を書いたので、そちらを参考(Phi2とPsi2が、ヘリシティ+1/-1の最低ウェイト状態)。空間回転とLorentz boostの作用には、$so(3)$の3次元表現に対応する項が加わるが、この表現の形を見る(AsirコードのT1,T2,T3を参照)と、twistor理論で扱われている、標準的な基底は、あんまりよくないという気がする(基底の取り方は、任意であって、良い悪いは完全に、人間の主観であるが、Riemann?Silberstein vectorのようなものを使って書く方がよいのかもしれない)。


あと、Minkowski空間上の場を考えると、共形変換群の作用は、純幾何学的に書けるので、Perelomovの意味でのgeneralized coherent states(lowest weight vectorへ$SU(2,2)$を作用させて得られる軌道の元)の具体的表示を得るのが容易になる(純表現論的に計算するのは、しんどい)。最低ウェイト条件から、特殊共形変換は、最低ウェイトベクトルを不変に保つし、並進は、純幾何学的な並進。

#このgeneralized coherent statesは、普通の量子光学で出てくるコヒーレント状態とは、全く別物で無関係。通常のコヒーレント状態は、様々な光子数状態の重ね合わせになっているけども、SU(2,2)-generalized coherent statesの方は、光子数は一個の状態(と考えるべきなのだと思う)。マクロなMaxwell方程式と考えてしまうと、物質がないのに電磁場が存在するというのは、奇妙な気がするし、その場合エネルギーや運動量は場の関数であるので、勝手に規格化することも物理的に変である



得られた最低ウェイト状態は、式としては、twistor理論の説明をしている文献で、時々見かけるものではある。それでも、こんな形で、ladder representationが実現されているというのは、ちょっと意外だった

/**** common stuffs ****/
load("noro_pd.rr");
load("noro_matrix.rr");

def assert(V){
   if(!V){error(V);}
   return 1;
}


def appdiff_rat(P , F , V){
   N = length(V)/2;
   if(type(P)==9){  /* 分散多項式表現の時 */
      DP = dp_dtop(P, V);
   }else if(type(P)<=2){
      DP = P;
   }else{
      return;
   }
   Terms = p_terms(DP , V , 2);
   for(Ret = 0 , Terms = Terms; Terms!=[] ; Terms=cdr(Terms)){
      T = car(Terms);
      if(type(F)==9){
           Tmp = dp_dtop(F , V);
      }else{
           Tmp = F;
      }
      if(T==1){
         Ret += p_nf(DP , V , V ,0)*Tmp;
      }else{
         C = p_nf(sdiv(DP-p_nf(DP,[T],V,0) , T) , V , V , 0);  /* 単項式Tの係数 */
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[N+I]);
            for(J = 0 ; J < Deg ; J++){
               Tmp = red(diff(Tmp , V[I]));
            }
         }
         for(I = 0 ; I < N ; I++){
            Deg = deg(T , V[I]);
            Tmp *= V[I]^Deg;
         }
         Ret += C*Tmp;
      }
   }
   return red(Ret);
}


def appdiff_complex(P , F , V){
    if(type(P)==9){
        DP = dp_dtop(P,V);
    }else{
        DP = P;
    }
    RP = real(DP);
    IP = imag(DP);
    if(type(F)<=2){
        RF = real(F);
        IF = imag(F);
    }else if(type(F)==3){ /*有理式*/
        DN_conj = conj(dn(F));
        if(DN_conj != dn(F)){
           TF = red(DN_conj*nm(F)/(DN_conj*dn(F)));
        }else{
           TF = F;
        }
        RF = red(real(nm(TF))/dn(TF));
        IF = red(imag(nm(TF))/dn(TF));
    }
    Ret = red( appdiff_rat(RP , RF , V) - appdiff_rat(IP,IF,V) );
    Ret += @i*red( appdiff_rat(IP,RF,V) + appdiff_rat(RP,IF,V) );
    return Ret;
}

/* 
微分作用Pを式Fに作用させる 
Vは変数のリスト(dp_dtopなどで使われるのと同様の形式)
*/
def appdiff(P,F,V){
    if(type(P)!=6){
        if(type(F)<=3){
           return appdiff_complex(P,F,V);
        }else if(type(F)==5){   /* vector */
           Ret = newvect(length(F));
           for(I=0;I<length(F);I++){
              Ret[I] = appdiff_complex(P, F[I] , V);
           }
           return Ret;
        }
    }else{ /*行列*/
        Sz = size(P);
        N = Sz[0];
        M = Sz[1];
        if(type(F)==5){
           Ret = newvect(N);
           for (I = 0 ; I < N ; I++){
              for(J = 0 ; J < M ; J++){
                Ret[I] += appdiff_complex(P[I][J] , F[J] , V);
              }
              Ret[I] = red(real(nm(Ret[I]))/dn(Ret[I])) + @i*red(imag(nm(Ret[I]))/dn(Ret[I]));
           }
           return Ret;
        }else if(type(F)==6){
            Sz2 = size(F);
            N2 = Sz2[0];
            M2 = Sz2[1];
            Ret = newmat(N,M2);
            for (I = 0 ; I < N ; I++){
                for(J = 0 ; J < M2 ; J++){
                    for(K = 0 ; K < M ; K++){
                       Ret[I][J] += appdiff_complex(P[I][K] , F[K][J] , V);
                    }
                }
            }
            return Ret;
        }
    }
}


def dp_weyl_matmul(A,B){
    Sz1 = size(A);
    Sz2 = size(B);
    L = Sz1[0];
    M = Sz1[1];
    N = Sz2[0];
    Ret = newmat(L,N);
    for(I = 0; I < L ; I++){
       for(J = 0 ; J < N ; J++){
          for(K = 0 ; K < M ; K++){
               Ret[I][J] += dp_weyl_mul(A[I][K] , B[K][J]);
          }
       }
    }
    return Ret;
}


V=[x0,x1,x2,x3,d0,d1,d2,d3];
X0=dp_ptod(x0,V);
X1=dp_ptod(x1,V);
X2=dp_ptod(x2,V);
X3=dp_ptod(x3,V);
D0=dp_ptod(d0,V);
D1=dp_ptod(d1,V);
D2=dp_ptod(d2,V);
D3=dp_ptod(d3,V);
C2 = D0*D0 - D1*D1 - D2*D2 - D3*D3;  /* Klein-Gordon operator */

/**** The end of common stuffs ****/


/**** scalar field stuffs ****/
Phi0 = 1/(x0^2-x1^2-x2^2-x3^2);    /* lowest weight state of helicity-0 massless particle */

/* infinitesimal conformal transformation for scalar field */
D = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3) + dp_ptod(1,V);
P0=D0;
P1=D1;
P2=D2;
P3=D3;
K0 = 2*dp_weyl_mul(X0,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0);  /* special conformal transformation */
K1 = -2*dp_weyl_mul(X1,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1); /* special conformal transformation */
K2 = -2*dp_weyl_mul(X2,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2); /* special conformal transformation */
K3 = -2*dp_weyl_mul(X3,D) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3); /* special conformal transformation */
M1 = dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2);
M2 = dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3);
M3 = dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1);
L1 = dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0);  /* Lorentz boost */
L2 = dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0);  /* Lorentz boost */
L3 = dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0);  /* Lorentz boost */


/* Chevalley generators for sl(4) and so(6) */
E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;

/* and other Cartan-Weyl Basis */
E4 = dp_weyl_mul(E1,E2) - dp_weyl_mul(E2,E1);   /* @i*(P1 + @i*P2)/2 */
E5 = dp_weyl_mul(E2,E3) - dp_weyl_mul(E3,E2);   /* -@i*(P1 - @i*P2)/2 */
E6 = dp_weyl_mul(E1,E5) - dp_weyl_mul(E5,E1);   /* (-P0+P3)/2 */
F4 = dp_weyl_mul(F1,F2) - dp_weyl_mul(F2,F1);   /* @i*(K1 - @i*K2)/2 */
F5 = dp_weyl_mul(F2,F3) - dp_weyl_mul(F3,F2);   /* -@i*(K1+@i*K2)/2 */
F6 = dp_weyl_mul(F1,F5) - dp_weyl_mul(F5,F1);   /* (K0+K3)/2 */
H4 = H1+H2;
H5 = H2+H3;
H6 = H1+H2+H3;

assert(appdiff(C2 , Phi0 , V)==0); /* Klein-Gordon equation */
assert(appdiff(F1,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(F2,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(F3,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(H1,Phi0,V)==0);     /* lowest weight condition */
assert(appdiff(H2,Phi0,V)==Phi0);  /* lowest weight condition */
assert(appdiff(H3,Phi0,V)==0);      /* lowest weight condition */
assert(appdiff(E1,Phi0,V)==0);     /* first singular vector */
assert(appdiff(E3,Phi0,V)==0);     /* second singular vector */


/**** The end of scalar field stuffs ****/


/**** "spin 1/2 fields" ****/
/* Pauli matrix */
S0 = newmat(2,2,[[1,0],[0,1]]);
S1 = newmat(2,2,[[0,1],[1,0]]);
S2 = newmat(2,2,[[0,-@i],[@i,0]]);
S3 = newmat(2,2,[[1,0],[0,-1]]);

/**** helicity +1/2 ****/
Phi = newvect(2,[(x1-@i*x2)/(x0^2-x1^2-x2^2-x3^2)^2 , (x0-x3)/(x0^2-x1^2-x2^2-x3^2)^2]);

DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);
D = (DS + dp_ptod(3/2,V))*S0;
P0 = D0*S0;
P1 = D1*S0;
P2 = D2*S0;
P3 = D3*S0;
M1 = S0*(dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2)) + dp_ptod(@i/2,V)*S1;
M2 = S0*(dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3)) + dp_ptod(@i/2,V)*S2;
M3 = S0*(dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1)) + dp_ptod(@i/2,V)*S3;
L1 = S0*(dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0)) + dp_ptod(-1/2,V)*S1;
L2 = S0*(dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0)) + dp_ptod(-1/2,V)*S2;
L3 = S0*(dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0)) + dp_ptod(-1/2,V)*S3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*S0 + 3*X0*S0 - X1*S1 - X2*S2 - X3*S3; 
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*S0 - 3*X1*S0 + X0*S1 - @i*X2*S3 + @i*X3*S2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*S0 - 3*X2*S0 + X0*S2 + @i*X1*S3 - @i*X3*S1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*S0 - 3*X3*S0 - @i*X1*S2 + @i*X2*S1 + X0*S3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3,Phi,V)==0*Phi );   /* Weyl equation */
assert( appdiff(F1,Phi,V)==0*Phi );
assert( appdiff(F2,Phi,V)==0*Phi );
assert( appdiff(F3,Phi,V)==0*Phi );
assert( appdiff(H1,Phi,V)+Phi==0*Phi );          /* j1==1/2 */
assert( appdiff(H2,Phi,V)-2*Phi==0*Phi );        /* d+j1+j2==2 */
assert( appdiff(H3,Phi,V)==0*Phi );              /* j2==0 */
/* singular vector vanishing */
assert( appdiff(E1,appdiff(E1,Phi,V),V)==0*Phi );
assert( appdiff(E3,Phi,V)==0*Phi ); 
assert( appdiff(E4,Phi,V)-appdiff(E2,appdiff(E1,Phi,V),V)==0*Phi );
assert( appdiff(E6 , appdiff(E2,Phi,V) , V) - appdiff(E4 , appdiff(E5,Phi,V) , V)==0*Phi);


/**** helicity -1/2 ****/
Psi = newvect(2, [-x0+x3 , x1+@i*x2])/(x0^2-x1^2-x2^2-x3^2)^2;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);
D = (DS + dp_ptod(3/2,V))*S0;
P0 = D0*S0;
P1 = D1*S0;
P2 = D2*S0;
P3 = D3*S0;
M1 = S0*(dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2)) + dp_ptod(@i/2,V)*S1;
M2 = S0*(dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3)) + dp_ptod(@i/2,V)*S2;
M3 = S0*(dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1)) + dp_ptod(@i/2,V)*S3;
L1 = S0*(dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0)) + dp_ptod(1/2,V)*S1;
L2 = S0*(dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0)) + dp_ptod(1/2,V)*S2;
L3 = S0*(dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0)) + dp_ptod(1/2,V)*S3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*S0 + 3*X0*S0 + X1*S1 + X2*S2 + X3*S3; 
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*S0 - 3*X1*S0 - X0*S1 - @i*X2*S3 + @i*X3*S2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*S0 - 3*X2*S0 - X0*S2 + @i*X1*S3 - @i*X3*S1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*S0 - 3*X3*S0 - @i*X1*S2 + @i*X2*S1 - X0*S3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , Psi , V)==0*Psi );   /* Weyl equation */
assert( appdiff(H1,Psi,V)==0*Psi );        /* j1==0 */
assert( appdiff(H2,Psi,V)-2*Psi==0*Psi );  /* d+j1+j2==2 */
assert( appdiff(H3,Psi,V)+Psi==0*Psi );    /* j2==1/2 */
assert( appdiff(F1,Psi,V)==0*Psi );
assert( appdiff(F2,Psi,V)==0*Psi );
assert( appdiff(F3,Psi,V)==0*Psi );
assert( appdiff(E1,Psi,V)==0*Psi );
assert( appdiff(E3,appdiff(E3,Psi,V),V)==0*Psi );
assert( appdiff(E5,Psi,V)+appdiff(E2,appdiff(E3,Psi,V),V)==0*Psi );
assert( appdiff(E6 , appdiff(E2,Psi,V) , V) - appdiff(E4 , appdiff(E5,Psi,V) , V)==0*Psi);


/**** "spin 1" fields ****/
T0 = dp_ptod(1,V)*newmat(3,3,[[1,0,0],[0,1,0],[0,0,1]]);
T1 = dp_ptod(1,V)*@i*newmat(3,3,[[0,-2,0],[-1,0,-1],[0,-2,0]])/2;
T2 = dp_ptod(1,V)*@i*newmat(3,3,[[0,2*@i,0],[-@i,0,@i],[0,-2*@i,0]])/2;
T3 = dp_ptod(1,V)*@i*newmat(3,3,[[-2,0,0],[0,0,0],[0,0,2]])/2;

/* (T1,T2,T3) is a representaion of so(3) */
assert(T1*T2-T2*T1==T3);
assert(T2*T3-T3*T2==T1);
assert(T3*T1-T1*T3==T2);


/**** helicity +1 ****/
Phi2 = newvect(3,[(x1-@i*x2)^2 , (x0-x3)*(x1-@i*x2) , (x0-x3)^2])/(x0^2-x1^2-x2^2-x3^2)^3;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);

D = (DS + 2)*T0;
P0 = D0*T0;
P1 = D1*T0;
P2 = D2*T0;
P3 = D3*T0;
L1 = (dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0))*T0 - @i*T1;
L2 = (dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0))*T0 - @i*T2;
L3 = (dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0))*T0 - @i*T3;
M1 = (dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2))*T0 - T1;
M2 = (dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3))*T0 - T2;
M3 = (dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1))*T0 - T3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*T0 + 4*X0*T0 - 2*@i*(X1*T1+X2*T2+X3*T3);
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*T0 - 4*X1*T0 + 2*@i*X0*T1 + 2*X2*T3 - 2*X3*T2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*T0 - 4*X2*T0 + 2*@i*X0*T2 - 2*X1*T3 + 2*X3*T1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*T0 - 4*X3*T0 + 2*X1*T2 - 2*X2*T1 + 2*@i*X0*T3;

E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );


/* SD-Maxwell equation */
assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Phi2[0],Phi2[1]]) , V)==newvect(2,[0,0]) );
assert( appdiff(D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Phi2[1],Phi2[2]]) , V)==newvect(2,[0,0]) );

assert( appdiff(H1,Phi2,V)+2*Phi2==0*Phi2 );
assert( appdiff(H2,Phi2,V)-3*Phi2==0*Phi2 );
assert( appdiff(H3,Phi2,V)==0*Phi2 );
assert( appdiff(F1,Phi2,V)==0*Phi2 );
assert( appdiff(F2,Phi2,V)==0*Phi2 );
assert( appdiff(F3,Phi2,V)==0*Phi2 );
assert( appdiff(E1,appdiff(E1,appdiff(E1,Phi2,V),V),V)==0*Phi2 );
assert( appdiff(E3,Phi2,V)==0*Phi2 );
assert( 2*appdiff(E4,Phi2,V)-appdiff(E2,appdiff(E1,Phi2,V),V)==0*Phi2 );
assert( appdiff(E6 , appdiff(E2,Phi2,V) , V) - appdiff(E4 , appdiff(E5,Phi2,V) , V)==0*Phi2);


/**** helicity -1 ****/
Psi2 = newvect(3,[(-x0+x3)^2 , (x1+@i*x2)*(-x0+x3) , (x1+@i*x2)^2])/(x0^2-x1^2-x2^2-x3^2)^3;
DS = dp_weyl_mul(X0,D0) + dp_weyl_mul(X1,D1) + dp_weyl_mul(X2,D2) + dp_weyl_mul(X3,D3);

D = (DS + 2)*T0;
P0 = D0*T0;
P1 = D1*T0;
P2 = D2*T0;
P3 = D3*T0;
L1 = (dp_weyl_mul(X0,D1) + dp_weyl_mul(X1,D0))*T0 + @i*T1;
L2 = (dp_weyl_mul(X0,D2) + dp_weyl_mul(X2,D0))*T0 + @i*T2;
L3 = (dp_weyl_mul(X0,D3) + dp_weyl_mul(X3,D0))*T0 + @i*T3;
M1 = (dp_weyl_mul(X2,D3) - dp_weyl_mul(X3,D2))*T0 - T1;
M2 = (dp_weyl_mul(X3,D1) - dp_weyl_mul(X1,D3))*T0 - T2;
M3 = (dp_weyl_mul(X1,D2) - dp_weyl_mul(X2,D1))*T0 - T3;
K0 = (2*dp_weyl_mul(X0,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D0))*T0 + 4*X0*T0 + 2*@i*(X1*T1+X2*T2+X3*T3);
K1 = (-2*dp_weyl_mul(X1,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D1))*T0 - 4*X1*T0 - 2*@i*X0*T1 + 2*X2*T3 - 2*X3*T2;
K2 = (-2*dp_weyl_mul(X2,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D2))*T0 - 4*X2*T0 - 2*@i*X0*T2 - 2*X1*T3 + 2*X3*T1;
K3 = (-2*dp_weyl_mul(X3,DS) - dp_weyl_mul(X0*X0-X1*X1-X2*X2-X3*X3,D3))*T0 - 4*X3*T0 + 2*X1*T2 - 2*X2*T1 - 2*@i*X0*T3;


E1 = (M1 + @i*M2)/2 + @i*(L1 + @i*L2)/2;
F1 = -(M1 - @i*M2)/2 - @i*(L1 - @i*L2)/2;
H1 = L3 - @i*M3;
E2 = (-P0 - P3)/2;
F2 = (K0 - K3)/2;
H2 = -(D + L3);
E3 = -(M1 - @i*M2)/2 + @i*(L1 - @i*L2)/2;
F3 = (M1 + @i*M2)/2 - @i*(L1 + @i*L2)/2;
H3 = L3 + @i*M3;
E4 = dp_weyl_matmul(E1,E2) - dp_weyl_matmul(E2,E1);
E5 = dp_weyl_matmul(E2,E3) - dp_weyl_matmul(E3,E2);
E6 = dp_weyl_matmul(E1,E5) - dp_weyl_matmul(E5,E1);

/* parts of sl(4) relation */
assert( dp_weyl_matmul(E1,F1) - dp_weyl_matmul(F1,E1)==H1 );
assert( dp_weyl_matmul(E2,F2) - dp_weyl_matmul(F2,E2)==H2 );
assert( dp_weyl_matmul(E3,F3) - dp_weyl_matmul(F3,E3)==H3 );
assert( dp_weyl_matmul(E2,F1)==dp_weyl_matmul(F1,E2) );
assert( dp_weyl_matmul(H1,E1) - dp_weyl_matmul(E1,H1)==2*E1 );
assert( dp_weyl_matmul(H2,E2) - dp_weyl_matmul(E2,H2)==2*E2 );

/* ASD Maxwell equation */
assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Psi2[0],Psi2[1]]) , V)==newvect(2,[0,0]) );
assert( appdiff(-D0*S0+D1*S1+D2*S2+D3*S3 , newvect(2,[Psi2[1],Psi2[2]]) , V)==newvect(2,[0,0]) );

assert( appdiff(H1,Psi2,V)==0*Psi2 );
assert( appdiff(H2,Psi2,V)-3*Psi2==0*Psi2 );
assert( appdiff(H3,Psi2,V)+2*Psi2==0*Psi2 );
assert( appdiff(F1,Psi2,V)==0*Psi2 );
assert( appdiff(F2,Psi2,V)==0*Psi2 );
assert( appdiff(F3,Psi2,V)==0*Psi2 );
assert( appdiff(E1,Psi2,V)==0*Psi2 );
assert( appdiff(E3,appdiff(E3,appdiff(E3,Psi2,V),V),V)==0*Psi2 );
assert( 2*appdiff(E5,Psi2,V)+appdiff(E2,appdiff(E3,Psi2,V),V)==0*Psi2 );
assert( appdiff(E6 , appdiff(E2,Psi2,V) , V) - appdiff(E4 , appdiff(E5,Psi2,V) , V)==0*Psi2);