FORM 4.0 (Jun 30 2012) 64-bits                   Run: Mon Feb  4 11:47:30 2013
    *   Title:  POINT CHARGE:
    *           energy-momentum tensor as a tempered distribution
    *    File:  energy-tensor.frm
    *   Tests:  energy-tensor-test.frm
    *  Author:  David N. Williams
    * License:  Creative Commons Attribution-Share Alike
    * Started:  August 28, 2012
    * Revised:  August 28, 2012
    *           September 3-5,7-11,28, 2012
    *           October 9,14, 2012
    *
    * All parts of this program not in the public domain are:
    *
    *   Copyright (C) 1987-1990, 2012 David N. Williams
    *
    * This work is licensed under the Creative Commons Attribution-
    * Share Alike 2.5 License.  To view a copy of this license, visit
    * http://creativecommons.org/licenses/by-sa/2.5/ or send a letter
    * to Creative Commons, 543 Howard Street, 5th Floor, San
    * Francisco, California, 94105, USA.
    *
    *   Input:  derivs.sav, derivs-ra.h
    *  Output:  energy-momentum.sav
    
    #if 0
    This program verifies the derivatives portion of a computation
    we finished before October 13, 1974, which showed how to realize
    the energy-momentum tensor for the retarded or advanced
    (ret/adv) electromagnetic field of a point charge as a tempered
    distribution.  We delayed publication, and ended up being
    scooped by the independent work of E. G. Peter Rowe, Phys. Rev.
    D18, 3639-54 (1978).  So we never submitted our result.
    
    Many of the precalculated derivative expressions in derivs.sav
    are not used.
    
    We use Gaussian units.  The source point and all vectors and
    tensors are retarded or advanced.  The Lorentz metric is (+---).
    Indices have the values 0,1,2,3.
    
    The source point and all vectors and tensors in the following
    are retarded or advanced.
    
    NOTATION
      c		 speed of light
      A		 Lorentz gauge four-vector Liénard-Wiechert potential
      F		 electromagnetic field tensor
      r		 light-cone four-displacement, from source to field point
      u		 four-velocity
      a		 four-acceleration
      j		 four-jerk
      s		 sgn(r⋅u) = sgn(r⁰)
      q		 charge
      f		 function of r⋅u
      ~t	 transpose of a second-rank tensor t
    
    Bracketed quantities [r], [u], [a], [f], etc., are used on the
    right-hand side of expressions that are to be differentiated by
    the "partial" procedure in derivs-ra.h.
    
    See derivs.frm for precalculated derivatives, such as ddf(mu,nu).
    
    All derivatives are understood to be valid when r and r⋅u are
    nonzero, i.e., off the world line.  Note that s = ±1 for
    ret/adv fields.  A factor qs is left out of A and its
    derivative.  Since it propagates through the derivative of A,
    occurrences of the field tensor F should be multiplied by it.
    That is, energy-momentum tensor results should be multiplied by
    (qs)² = q².
    #endif
    
    *** DECLARATIONS
    
    #include derivs-ra.h
    #-
    
    off statistics;
    I mu,nu,la,ro;
    * tensors for pretty printing
    CF [r,u],[r,a],[rr],[{r,u}],[{r,a}],[uu],[{u,a}],[aa];
    
    Load derivs.sav;
 dr loaded
 du loaded
 da loaded
 dj loaded
 [d(r.u)] loaded
 [d(r.a)] loaded
 [d.r] loaded
 [d.u] loaded
 [d.a] loaded
 [d.j] loaded
 [r.dr] loaded
 [r.du] loaded
 [r.da] loaded
 [r.dj] loaded
 [u.dr] loaded
 [u.du] loaded
 [u.da] loaded
 [a.dr] loaded
 [a.du] loaded
 [a.da] loaded
 df loaded
 df1 loaded
 [r.d(r.u)] loaded
 [r.d(r.a)] loaded
 [u.d(r.u)] loaded
 [u.d(r.a)] loaded
 [a.d(r.u)] loaded
 [a.d(r.a)] loaded
 [r.df] loaded
 [u.df] loaded
 dfr loaded
 dfu loaded
 dfa loaded
 dfj loaded
 df1r loaded
 df1u loaded
 df1a loaded
 df1j loaded
 drr loaded
 dru loaded
 dra loaded
 drj loaded
 duu loaded
 dua loaded
 duj loaded
 daa loaded
 daj loaded
 dfrr loaded
 dfru loaded
 dfra loaded
 dfrj loaded
 dfuu loaded
 dfua loaded
 dfuj loaded
 dfaa loaded
 dfaj loaded
 [d.fr] loaded
 [d.fu] loaded
 [d.fa] loaded
 [d.fj] loaded
 [d.frr] loaded
 [d.fru] loaded
 [d.fra] loaded
 [d.fur] loaded
 [d.fuu] loaded
 [d.fua] loaded
 [d.far] loaded
 [d.fau] loaded
 [d.faa] loaded
 ddf loaded
 ddr loaded
 ddu loaded
 dda loaded
 ddfr loaded
 ddfu loaded
 ddfa loaded
 ddrr loaded
 ddru loaded
 ddra loaded
 dduu loaded
 ddua loaded
 ddaa loaded
 ddfrr loaded
 ddfru loaded
 ddfra loaded
 ddfuu loaded
 ddfua loaded
 ddfaa loaded
 dddfr loaded
    .global
    
    *** PROCEDURES
    
    * debugging
    #procedure try
    print +s;
    .end
    #endprocedure
    
    * NOTE:  no dummies in these two procedures
    #procedure idantisym
      Id r(mu)*u(nu) = [r,u](mu,nu) + r(nu)*u(mu);
      Id r(mu)*a(nu) = [r,a](mu,nu) + r(nu)*a(mu);
      Id u(mu)*a(nu) = [u,a](mu,nu) + u(nu)*a(mu);
    #endprocedure
    *
    #procedure idsym
      Id r(mu)*r(nu) = [rr](mu,nu);
      Id r(mu)*u(nu) = [{r,u}](mu,nu) - r(nu)*u(mu);
      Id r(mu)*a(nu) = [{r,a}](mu,nu) - r(nu)*a(mu);
      Id u(mu)*u(nu) = [uu](mu,nu);
      Id a(mu)*a(nu) = [aa](mu,nu);
    #endprocedure
    
    
    *** MASTER SHEET
    
    *
    * The tensors in this section are from the the 1974 master
    * sheet.  Unfortunately it uses notation like (A⋅A)² for
    * c⁴/(r⋅u)⁴, which is misleading, because the energy-momentum
    * tensor is actually quadratic in A.  We keep that notation in
    * the commented formulas to facilitate comparison with the
    * master sheet.
    *
    
    * A = u/(r⋅u) (omits qs)
    G A(mu) = f*u(mu);
    #call idf
    #call idmetric
    print +s;
    .store

   A(mu) =
       + u(mu)*r.u^-1
      ;

    
    * F = ∂A - ~∂A
    G F(mu,nu) = dfu(mu,nu) - dfu(nu,mu);
    #call idf
    #call idmetric
    *#call idantisym
    *B [r,u],[r,a];
    print +s;
    .store

   F(mu,nu) =
       - r(mu)*u(nu)*r.u^-3*r.a
       + r(mu)*u(nu)*r.u^-3*c^2
       + r(mu)*a(nu)*r.u^-2
       + r(nu)*u(mu)*r.u^-3*r.a
       - r(nu)*u(mu)*r.u^-3*c^2
       - r(nu)*a(mu)*r.u^-2
      ;

    
    * t₁ = 1/2 ∂⋅∂(A⋅A) g = 1/2 c²∂⋅∂[1/(r⋅u)²)] g
    G t1(mu,nu) = 1/2*c^2*ddf(la,la)*g(mu,nu);
    #call idfp2
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t1(mu,nu) =

       + g(mu,nu) * (
          + 2*r.u^-4*r.a*c^2
          - r.u^-4*c^4
          );

    
    * t₂ = ∂∂(A⋅A) = c²∂∂[1/(r⋅u)²)]
    G t2(mu,nu) = c^2*ddf(mu,nu);
    #call idfp2
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t2(mu,nu) =

       + g(mu,nu) * (
          - 2*r.u^-4*r.a*c^2
          + 2*r.u^-4*c^4
          )

       + [rr](mu,nu) * (
          - 16*r.u^-6*r.a*c^4
          + 8*r.u^-6*r.a^2*c^2
          + 8*r.u^-6*c^6
          - 2*r.u^-5*r.j*c^2
          )

       + [{r,u}](mu,nu) * (
          + 8*r.u^-5*r.a*c^2
          - 8*r.u^-5*c^4
          )

       + [{r,a}](mu,nu) * (
          - 2*r.u^-4*c^2
          )

       + 6*[uu](mu,nu)*r.u^-4*c^2
         ;

    
    * t₃ =3/4 ∂⋅∂[(A⋅A)²rr] = 3/4 c⁴∂⋅∂[rr/(r⋅u)⁴]
    G t3(mu,nu) = 3/4*c^4*ddfrr(la,la,mu,nu);
    #call idfp4
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t3(mu,nu) =

       + g(mu,nu) * (
          + 3/2*r.u^-4*c^4
          )

       + [rr](mu,nu) * (
          + 6*r.u^-6*r.a*c^4
          + 3*r.u^-6*c^6
          )

       + [{r,u}](mu,nu) * (
          - 3*r.u^-5*c^4
          );

    
    * t₄ = -3/4 ∂⋅∂[A⋅A(rA + Ar)] = -3/4 c²∂⋅∂[(ru + ur)/(r⋅u)³]
     
    G t4(mu,nu) = - 3/4*c^2*( ddfru(la,la,mu,nu) + ddfru(la,la,nu,mu) );
    #call idfp3
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t4(mu,nu) =

       + [{r,u}](mu,nu) * (
          - 9/2*r.u^-5*r.a*c^2
          )

       + [{r,a}](mu,nu) * (
          + 3/2*r.u^-4*c^2
          )

       + 3*[uu](mu,nu)*r.u^-4*c^2
         ;

    
    * t₅ = -3/4 ∂⋅∂[∂(A⋅Ar) + ~∂(A⋅Ar)]
    *    = -3/4 c²∂⋅∂{∂[r/(r⋅u)²] + ~∂[r/(r⋅u)²]}
    G t5(mu,nu) = -3/4*c^2*( dddfr(la,la,mu,nu) + dddfr(la,la,nu,mu) );
    #call idfp2
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t5(mu,nu) =

       + g(mu,nu) * (
          - 3*r.u^-4*c^4
          )

       + [rr](mu,nu) * (
          + 12*r.u^-6*r.a*c^4
          - 12*r.u^-6*c^6
          )

       + [{r,u}](mu,nu) * (
          - 9/2*r.u^-5*r.a*c^2
          + 12*r.u^-5*c^4
          )

       + [{r,a}](mu,nu) * (
          + 3/2*r.u^-4*c^2
          )

       - 9*[uu](mu,nu)*r.u^-4*c^2
         ;

    
    * w₅ = -3/4 [∂(A⋅Ar) + ~∂(A⋅Ar)]
    *    = -3/4 c²{∂[r/(r⋅u)²] + ~∂[r/(r⋅u)²]}
    * ∂⋅∂ w₅ = t₅
    G w5(mu,nu) = -3/4*c^2*( dfr(mu,nu) + dfr(nu,mu) );
    #call idfp2
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   w5(mu,nu) =

       + g(mu,nu) * (
          - 3/2*r.u^-2*c^2
          )

       + [rr](mu,nu) * (
          + 3*r.u^-4*r.a*c^2
          - 3*r.u^-4*c^4
          )

       + [{r,u}](mu,nu) * (
          + 9/4*r.u^-3*c^2
          );

    
    * This does not vanish.
    G [d.w5](mu) = -3/4*c^2*( ddfr(nu,nu,mu) + ddfr(nu,mu,nu) );
    #call idfp2
    #call idmetric
    print +s;
    .store

   [d.w5](mu) =
       + 3/2*r(mu)*r.u^-4*r.a*c^2
       - 3*r(mu)*r.u^-4*c^4
       + 3*u(mu)*r.u^-3*c^2
      ;

    
    * t₆ = rr[-a⋅a - 9c²(r⋅a)²/(r⋅u)² + 2c²r⋅j/r⋅u]/(r⋅u)⁴
    G t6(mu,nu) = r(mu)*r(nu)*
      ( -a.a - 9*c^2*r.a^2/r.u^2 + 2*c^2*r.j/r.u )/r.u^4;
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t6(mu,nu) =

       + [rr](mu,nu) * (
          - 9*r.u^-6*r.a^2*c^2
          + 2*r.u^-5*r.j*c^2
          - r.u^-4*a.a
          );

    
    * Maxwell energy-momentum tensor, absent a factor 1/4π.
    * The extra substitution for FF is a temporary workaround for a
    * bug in FORM, suggested by Jos Vermaseren.
    CF FF;
    G TCEM(mu,nu) = FF(mu,la)*FF(la,nu) + 1/4*g(mu,nu)*FF(la,ro)*FF(la,ro);
    Id FF(mu?,nu?) = F(mu,nu);
    #call idmetric
    #call idsym
    print +s;
    .store

   TCEM(mu,nu) =
       - 1/2*g(mu,nu)*r.u^-4*c^4
       + 2*[rr](mu,nu)*r.u^-6*r.a*c^4
       - [rr](mu,nu)*r.u^-6*r.a^2*c^2
       - [rr](mu,nu)*r.u^-6*c^6
       - [rr](mu,nu)*r.u^-4*a.a
       - [{r,u}](mu,nu)*r.u^-5*r.a*c^2
       + [{r,u}](mu,nu)*r.u^-5*c^4
       + [{r,a}](mu,nu)*r.u^-4*c^2
      ;

    
    
    *** REARRANGEMENT
    
    *
    * The tensors in this section are from the master sheet
    * supplementary page.  They rearrange the energy-momentum tensor
    * into a sum of terms with vanishing four-divergence off the
    * world line, which are used in our hand calculation of the
    * divergence on the world line.  The Tensor wi, wii,and wiii are
    * factors whose divergences vanish before a d'Alembertian is
    * applied.
    *
    
    * wi = -9/4 rr (A⋅A)²        ti = ⃞ wi
    G [d.wi](mu) = -9/4*c^4*dfrr(nu,nu,mu);
    #call idfp4
    #call idmetric
    print +s;
    .store

   [d.wi](mu) = 0;

    
    * wii = 3/2 rA (A⋅A)        tii =  ⃞ wiii
    G [d.wii](mu) = 3/2*c^4*dfru(nu,nu,mu);
    #call idfp3
    #call idmetric
    print +s;
    .store

   [d.wii](mu) = 0;

    
    * waiii = 3/2 Ar (A⋅A)       taiii =  ⃞ waiii 
    G [d.waiii](mu) = 3/2*c^2*dfru(nu,mu,nu);
    #call idfp3
    #call idmetric
    print +s;
    .store

   [d.waiii](mu) =
       - 3*r(mu)*r.u^-4*r.a*c^2
      ;

    
    * wbiii = 3 rr (A⋅A)² r⋅a/c²   tbiii =  ⃞ wbiii
    G [d.wbiii](mu) = 3*[r](nu)*[r](mu)*[f]*[r](la)*[a](la)*c^2;
    #call partial(nu)
    #call idvec
    #call idmetric
    #call idfp4
    print +s;
    .store

   [d.wbiii](mu) =
       + 3*r(mu)*r.u^-4*r.a*c^2
      ;

    
    * wiii = waiii + wbiii      tiii =  ⃞ wiii
    G [d.wiii](mu) = [d.waiii](mu) + [d.wbiii](mu);
    print;
    .store

   [d.wiii](mu) = 0;

    
    * Note that all three pieces of t6 = t6a + t6b + t6c have
    * vanishing divergence off the world line.
    
    * t6a = -rr a⋅a/(r⋅u)⁴
    G [d.t6a](mu) = -[r](nu)*[r](mu)*[f]*[a](la)*[a](la);
    #call partial(nu)
    #call idvec
    #call idmetric
    #call idfp4
    print;
    .store

   [d.t6a](mu) = 0;

    
    * t6b =  -9 rr c²(r⋅a)²/(r⋅u)⁶
    G [d.t6b](mu)= -9*[r](nu)*[r](mu)*[f]*c^2*[r](la)*[a](la)*[r](ro)*[a](ro);
    #call partial(nu)
    #call idvec
    #call idmetric
    #call idfp6
    print;
    .store

   [d.t6b](mu) = 0;

    
    * t6c =   2 rr c² r⋅j/r⋅u⁵ 
    G [d.t6c](mu) = 2*[r](nu)*[r](mu)*[f]*c^2*[r](la)*[j](la);
    #call partial(nu)
    #call idvec
    #call idmetric
    #call idfp5
    print;
    .store

   [d.t6c](mu) = 0;

    
    * t₇ = -g ⃞(A⋅A) + ∂∂(A⋅A) = (-g ⃞ + ∂∂) c²/(r⋅u)² 
    * ∂⋅t₇ = 0 by inspection of the above line
    G t7(mu,nu) = - g(mu,nu)*ddfuu(la,la,ro,ro) + ddfuu(mu,nu,la,la);
    #call idfp2
    #call idmetric
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   t7(mu,nu) =

       + g(mu,nu) * (
          - 6*r.u^-4*r.a*c^2
          + 4*r.u^-4*c^4
          )

       + [rr](mu,nu) * (
          - 16*r.u^-6*r.a*c^4
          + 8*r.u^-6*r.a^2*c^2
          + 8*r.u^-6*c^6
          - 2*r.u^-5*r.j*c^2
          )

       + [{r,u}](mu,nu) * (
          + 8*r.u^-5*r.a*c^2
          - 8*r.u^-5*c^4
          )

       + [{r,a}](mu,nu) * (
          - 2*r.u^-4*c^2
          )

       + 6*[uu](mu,nu)*r.u^-4*c^2
         ;

    
    * tᵢ = -9/4 ⃞ [rr c⁴/(r⋅u)⁴] 
    G ti(mu,nu) = -9/4*[r](mu)*[r](nu)*[f]*c^4;
    #call partial(la)
    #call partial(la)
    #call idvec
    #call idmetric
    #call idfp4
    print +s;
    .store

   ti(mu,nu) =
       - 18*r(mu)*r(nu)*r.u^-6*r.a*c^4
       - 9*r(mu)*r(nu)*r.u^-6*c^6
       + 9*r(mu)*u(nu)*r.u^-5*c^4
       + 9*r(nu)*u(mu)*r.u^-5*c^4
       - 9/2*g(mu,nu)*r.u^-4*c^4
      ;

    
    * tᵢᵢ = 3/2 ⃞ [ru c²/(r⋅u)³]
    G tii(mu,nu) = 3/2*[r](mu)*[u](nu)*[f]*c^2;
    #call partial(la)
    #call partial(la)
    #call idvec
    #call idmetric
    #call idfp3
    print +s;
    .store

   tii(mu,nu) =
       + 9*r(mu)*u(nu)*r.u^-5*r.a*c^2
       - 3*r(mu)*a(nu)*r.u^-4*c^2
       - 3*u(mu)*u(nu)*r.u^-4*c^2
      ;

    
    * taᵢᵢᵢ = 3/2 ⃞ [ur c²/(r⋅u)³]
    G taiii(mu,nu) = 3/2*[u](mu)*[r](nu)*[f]*c^2;
    #call partial(la)
    #call partial(la)
    #call idvec
    #call idmetric
    #call idfp3
    print +s;
    .store

   taiii(mu,nu) =
       + 9*r(nu)*u(mu)*r.u^-5*r.a*c^2
       - 3*r(nu)*a(mu)*r.u^-4*c^2
       - 3*u(mu)*u(nu)*r.u^-4*c^2
      ;

    
    * tbᵢᵢᵢ = 3 ⃞ [rr c² r⋅a/(r⋅u)⁴]
    G tbiii(mu,nu) = 3*[r](mu)*[r](nu)*c^2*[f]*[r](ro)*[a](ro);
    #call partial(la)
    #call partial(la)
    #call idvec
    #call idmetric
    #call idfp4
    print +s;
    .store

   tbiii(mu,nu) =
       + 36*r(mu)*r(nu)*r.u^-6*r.a*c^4
       - 18*r(mu)*u(nu)*r.u^-5*r.a*c^2
       + 6*r(mu)*a(nu)*r.u^-4*c^2
       - 18*r(nu)*u(mu)*r.u^-5*r.a*c^2
       + 6*r(nu)*a(mu)*r.u^-4*c^2
       + 6*g(mu,nu)*r.u^-4*r.a*c^2
      ;

    
    G T(mu,nu) = t7(mu,nu) + ti(mu,nu)  + tii(mu,nu)
       + taiii(mu,nu)  + tbiii(mu,nu)  + t6(mu,nu) 
    ;
    #call idsym
    B g,[rr],[{r,u}],[{r,a}];
    print +s;
    .store

   T(mu,nu) =

       + g(mu,nu) * (
          - 1/2*r.u^-4*c^4
          )

       + [rr](mu,nu) * (
          + 2*r.u^-6*r.a*c^4
          - r.u^-6*r.a^2*c^2
          - r.u^-6*c^6
          - r.u^-4*a.a
          )

       + [{r,u}](mu,nu) * (
          - r.u^-5*r.a*c^2
          + r.u^-5*c^4
          )

       + [{r,a}](mu,nu) * (
          + r.u^-4*c^2
          );

    
    save energy-tensor.sav;
    .end
  0.02 sec out of 0.06 sec