$SPAD/src/input biquat.input

D. Cyganski December 3, 2016

Abstract

Contents

— * —

)set break resume
)spool biquat.output
)set message test on
)set message auto off
)clear all

———-

1 Hamiltonian Biquaternions

Define the type Q of Hamiltonian biquaternions (not to be confused with Clifford biquaternions which are an entirely different object)

— * —

--S 1 of 43
C:=Complex Expression Integer
--R
--R
--R   (1)  Complex(Expression(Integer))
--R                                                                 Type: Domain
--E 1
--S 2 of 43
Q:=Quaternion C
--R
--R
--R   (2)  Quaternion(Complex(Expression(Integer)))
--R                                                                 Type: Domain
--E 2
--S 3 of 43
q:Q:=quatern(q0,q1,q2,q3)
--R
--R
--R   (3)  q0 + q1 i + q2 j + q3 k
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 3

———-

Define a function that takes a complex list (parameter l) into a quaternion — * —

--S 4 of 43
qlist(l:List C):Q==quatern(1.1,1.2,1.3,1.4)
--R
--R   Function declaration qlist : List(Complex(Expression(Integer))) ->
--R      Quaternion(Complex(Expression(Integer))) has been added to
--R      workspace.
--R                                                                   Type: Void
--E 4

———-

Define a function that takes a quaternion into a list — * —

--S 5 of 43
listq(x:Q):List C == [real x, imagI x, imagJ x, imagK x]
--R
--R   Function declaration listq : Quaternion(Complex(Expression(Integer))
--R      ) -> List(Complex(Expression(Integer))) has been added to
--R      workspace.
--R                                                                   Type: Void
--E 5

———-

Define a function that takes a biquat into a matrix — * —

--S 6 of 43
matrixq(x:Q):Matrix C == matrix _
             [[real x + %i*imagI(x), imagJ x + %i*imagK(x)],_
             [-imagJ(x) + %i*imagK(x), real x - %i*imagI(x)]]
--R
--R   Function declaration matrixq : Quaternion(Complex(Expression(Integer
--R      ))) -> Matrix(Complex(Expression(Integer))) has been added to
--R      workspace.
--R                                                                   Type: Void
--E 6

———-

2 The Pauli basis Biquaternions

The quaternion package in axiom is based upon the Hamiltonian basis, i, j, k, in contrast to the Pauli basis used in the Morgan papers which make certain operations correspond more closely (in the sense of visualization as ordering of components is more natural so that classical operations can be easily identified) to vector calculus functions and connect with more directly with various spinor QM formulations. To both follow these papers and obtain these visualization benefits, we will define and use the Pauli basis also.

The Pauli basis in terms of the Hamiltonian quaternion basis — * —

--S 7 of 43
sig0:=quatern(1,0,0,0)::Q
--R
--R
--R   (7)  1
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 7
--S 8 of 43
sig1:=%i*quatern(0,0,0,1)::Q
--R
--R
--R   (8)  %i k
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 8
--S 9 of 43
sig2:=%i*quatern(0,0,1,0)::Q
--R
--R
--R   (9)  %i j
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 9
--S 10 of 43
sig3:=-%i*quatern(0,1,0,0)::Q
--R
--R
--R   (10)  - %i i
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 10

———-

For purposes of manipulation and display it would be best if we either

1. Copy the quaternion source code and edit it to use σ_{i} as symbols for the basis rather than i, j, k and change the definition of multiplication to correspond to this Pauli basis so that all current display and entry functions are flipped to the new basis

2. Use the σ basis above with the existing package but define a new display function, say sigrep, that can show us the results in Pauli form at will, but otherwise show it in Hamiltonian basis

While 2 is not as natural and convenient, it is easier to code and would automatically provide a development that could be seen in the visually appealing Pauli form, or in the standard Hamiltonian form of biquaternions.

After a quick review and finding no fast path to 2, and not wishing to get into the mechanics of compiling an axiom package instead of working on the physics, I have taken the yet easier path and simply defined the simplest display function that gives me access to seeing Pauli basis results: a list of the four Pauli basis coefficients.

Define a function that produces the Pauli basis representation of the biquaternion. — * —

--S 11 of 43
siglist(x:Q):List C == [real x, -imagK(x)*%i, -imagJ(x)*%i, %i*imagI(x)]
--R
--R   Function declaration siglist : Quaternion(Complex(Expression(Integer
--R      ))) -> List(Complex(Expression(Integer))) has been added to
--R      workspace.
--R                                                                   Type: Void
--E 11

———-

3 Biquaternion Calculus

Define the quaternion derivative (Morgan, 2001, Eq. 2)

— * —

--S 12 of 43
D(q:Q,x:Symbol,y:Symbol,z:Symbol):Q==sig1*D(q,x)+sig2*D(q,y)+sig3*D(q,z)
--R
--R   Function declaration D : (Quaternion(Complex(Expression(Integer))),
--R      Symbol,Symbol,Symbol) -> Quaternion(Complex(Expression(Integer)))
--R      has been added to workspace.
--R                                                                   Type: Void
--E 12

———-

For testing the derivative we define this set of operators — * —

--S 13 of 43
Ft:=operator 'Ft
--R
--R
--R   (13)  Ft

--R                                                          Type: BasicOperator
--E 13
--S 14 of 43
Fx:=operator 'Fx
--R
--R
--R   (14)  Fx
--R                                                          Type: BasicOperator
--E 14
--S 15 of 43
Fy:=operator 'Fy
--R
--R
--R   (15)  Fy
--R                                                          Type: BasicOperator
--E 15
--S 16 of 43
Fz:=operator 'Fz
--R
--R
--R   (16)  Fz
--R                                                          Type: BasicOperator
--E 16

———-

Now form a general quaternion which is a function of x, y, z — * —

--S 17 of 43
F:Q:=Ft(x,y,z)*sig0+Fx(x,y,z)*sig1+Fy(x,y,z)*sig2+Fz(x,y,z)*sig3
--R
--R
--R   (17)  Ft(x,y,z) - Fz(x,y,z)%i i + Fy(x,y,z)%i j + Fx(x,y,z)%i k
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 17

———-

In the Pauli basis the derivative of this biquat should produce (Morgan 2001, eq 1):

  D(Ft+F.sigma)=div(F)+(grad(Ft)+%i*curl(F)).sigma

— * —

--S 18 of 43
siglist(D(F,x,y,z))
--R
--R   Compiling function D with type (Quaternion(Complex(Expression(
--R      Integer))),Symbol,Symbol,Symbol) -> Quaternion(Complex(Expression
--R      (Integer)))
--R   Compiling function siglist with type Quaternion(Complex(Expression(
--R      Integer))) -> List(Complex(Expression(Integer)))
--R
--R   (18)
--R   [Fz  (x,y,z) + Fy  (x,y,z) + Fx  (x,y,z),
--R      ,3            ,2            ,1
--R    Ft  (x,y,z) + (Fz  (x,y,z) - Fy  (x,y,z))%i,
--R      ,1             ,2            ,3
--R    Ft  (x,y,z) + (- Fz  (x,y,z) + Fx  (x,y,z))%i,
--R      ,2               ,1            ,3
--R    Ft  (x,y,z) + (Fy  (x,y,z) - Fx  (x,y,z))%i]
--R      ,3             ,1            ,2
--R                                     Type: List(Complex(Expression(Integer)))
--E 18

———-

4 Biquaternion Exponential

Define a function that evaluates the Biquat exponential by beginning with the biquat rotation.

First let’s introduce some refinement to the language we will use. In many texts and papers, the i, j, k or σ_{1}, σ_{2}, σ_{3} components of the biquaternion are called the imaginary components while the scalar or σ_{0} component are called real. This introduces vast confusion when dealing with biquaternions in which there is a “true” imaginary unit, %i being used. We will adopt a variation on the language used by Gsponer and others. A biquat [a, b, c, d] has scalar component a which may have itself an imaginary component. The component [0, b, c, d] is its vector component which again may be complex:

 [0,b,c,d]=[0,Re(b),Re(c),Re(d)]+[0,Im(b),Im(c),Im(d)]=Re(Vec(q))+Im(Vec(q))

Thus, for us, a “real quaternion” has real scalar and vector components and does not refer to a quaternion with only a scalar component as it would in the language adopted by the quatern package in Axiom.

The rotation is a function of real number, theta and a unit norm, real vector quaternion. It’s form is that of a simple complex quaternion exponential, R(θ,q)=e^{a} where a is an Imaginary Vector quaternions are like exponentials of pure imaginary numbers, resulting in each case in a unit norm object of the same type. For many of our purposes in the representation of relativistic and quantum mechanics it is the only type of quaternion for which we will have to evaluate the exponential.

Define a biquaternion rotation operator that takes a biquat through a rotation of θ radians about the axis defined by the unit q biquat (Morgan 2001, Eq 3).

— * —

--S 19 of 43
rot(theta:Expression Integer,q:Q):Q==cos(theta/2)-%i::Q*q*sin(theta/2)
--R
--R   Function declaration rot : (Expression(Integer),Quaternion(Complex(
--R      Expression(Integer)))) -> Quaternion(Complex(Expression(Integer))
--R      ) has been added to workspace.
--R                                                                   Type: Void
--E 19

———-

The rotation is a basis for defining the general exponential, since we can always extract the unit vector corresponding to a given biquaternion, all we need is the biquaternion abs operation, not provided by axiom (note, many texts and papers label the abs which is a complex valued norm^{2} by the name “norm” – very bad form – I’ve adopted abs as the name to avoid this confusion.

— * —

--S 20 of 43
((x:Q)/(y:Q)):Q == x*inv(y)
--R
--R   Function declaration ?/? : (Quaternion(Complex(Expression(Integer)))
--R      ,Quaternion(Complex(Expression(Integer)))) -> Quaternion(Complex(
--R      Expression(Integer))) has been added to workspace.
--R                                                                   Type: Void
--E 20
--S 21 of 43
abs(q:Q):C == sqrt((q*conjugate(q))::C)
--R
--R   Function declaration abs : Quaternion(Complex(Expression(Integer)))
--R       -> Complex(Expression(Integer)) has been added to workspace.
--R                                                                   Type: Void
--E 21
--S 22 of 43
exp(q:Q):Q == (_
  q-conjugate(q)=0 => exp( (q+conjugate(q))::C/2)$C * sig0; _
  exp( (q+conjugate(q))::C/2)$C * (sig0*cos(abs(q)) +_
  (q-conjugate(q))/abs(q-conjugate(q))*sin(abs(q))) )
--R
--R   Function declaration exp : Quaternion(Complex(Expression(Integer)))
--R       -> Quaternion(Complex(Expression(Integer))) has been added to
--R      workspace.

--R                                                                   Type: Void
--E 22

———-

If I’ve defined these correctly, then the rotation about the x axis defined by qx below by 2 radians should give the same answer as exponentiation of %i*qx (not a very complete test).

— * —

--S 23 of 43
qx:=sig1
--R
--R
--R   (23)  %i k
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 23
--S 24 of 43
mm:=siglist(rot(2,qx))
--R
--R   Compiling function / with type (Quaternion(Complex(Expression(
--R      Integer))),Quaternion(Complex(Expression(Integer)))) ->
--R      Quaternion(Complex(Expression(Integer)))
--R   There are 2 exposed and 6 unexposed library operations named cos
--R      having 1 argument(s) but none was determined to be applicable.
--R      Use HyperDoc Browse, or issue
--R                               )display op cos
--R      to learn more about the available operations. Perhaps
--R      package-calling the operation or using coercions on the arguments
--R      will allow you to apply the operation.
--R   Cannot find a definition or applicable library operation named cos
--R      with argument type(s)
--R                  Quaternion(Complex(Expression(Integer)))
--R
--R      Perhaps you should use "@" to indicate the required return type,
--R      or "$" to specify which version of the function you need.
--R   AXIOM will attempt to step through and interpret the code.
--R   Compiling function rot with type (Expression(Integer),Quaternion(
--R      Complex(Expression(Integer)))) -> Quaternion(Complex(Expression(
--R      Integer)))
--R
--R   (24)
--R                2         3                  3          2
--R    cos(1)sin(1)  + cos(1)  + cos(1)   sin(1)  + (cos(1)  - 1)sin(1)
--R   [-------------------------------- + ----------------------------- %i,
--R                  2          2                     2          2
--R           2sin(1)  + 2cos(1)               2sin(1)  + 2cos(1)
--R                  2         3                    3            2
--R    - cos(1)sin(1)  - cos(1)  + cos(1)   - sin(1)  + (- cos(1)  - 1)sin(1)

--R    ---------------------------------- + --------------------------------- %i,
--R                   2          2                        2          2
--R            2sin(1)  + 2cos(1)                  2sin(1)  + 2cos(1)
--R    0, 0]
--R                                     Type: List(Complex(Expression(Integer)))
--E 24
--S 25 of 43
nn:=siglist(exp(-%i::Q*qx))
--R
--R   There are 2 exposed and 7 unexposed library operations named exp
--R      having 1 argument(s) but none was determined to be applicable.
--R      Use HyperDoc Browse, or issue
--R                               )display op exp
--R      to learn more about the available operations. Perhaps
--R      package-calling the operation or using coercions on the arguments
--R      will allow you to apply the operation.
--R   Cannot find a definition or applicable library operation named exp
--R      with argument type(s)
--R                  Quaternion(Complex(Expression(Integer)))
--R
--R      Perhaps you should use "@" to indicate the required return type,
--R      or "$" to specify which version of the function you need.
--R   AXIOM will attempt to step through and interpret the code.
--R   Compiling function exp with type Quaternion(Complex(Expression(
--R      Integer))) -> Quaternion(Complex(Expression(Integer)))
--R   Compiling function abs with type Quaternion(Complex(Expression(
--R      Integer))) -> Complex(Expression(Integer))
--R
--R   (25)
--R                2         3                  3          2
--R    cos(1)sin(1)  + cos(1)  + cos(1)   sin(1)  + (cos(1)  - 1)sin(1)
--R   [-------------------------------- + ----------------------------- %i,
--R                  2          2                     2          2
--R           2sin(1)  + 2cos(1)               2sin(1)  + 2cos(1)
--R                  2         3                    3            2
--R    - cos(1)sin(1)  - cos(1)  + cos(1)   - sin(1)  + (- cos(1)  - 1)sin(1)
--R    ---------------------------------- + --------------------------------- %i,
--R                   2          2                        2          2
--R            2sin(1)  + 2cos(1)                  2sin(1)  + 2cos(1)
--R    0, 0]
--R                                     Type: List(Complex(Expression(Integer)))
--E 25
--S 26 of 43
(mm=nn)@Boolean
--R
--R
--R   (26)  true
--R                                                                Type: Boolean

--E 26

———-

To express a proof of equality such as:

  rot(theta,q) = exp((-theta/2)*%i*q)

for arbitrary real θ and biquaternion real using vector q as one would in Maple, we need to express the fact that q is a vector quantity with real coefficients in the Pauli basis representing a unit magnitude vector. One way to represent this is to create a quaternion called qnv, representing a general normalized (unit length) vector part only biquat with two independent variables q1 and q2 representing its only degrees of freedom.

— * —

--S 27 of 43
qnv:=q1*sig1+q2*sig2+sqrt(1-q1^2-q2^2)*sig3
--R
--R
--R            +---------------+
--R            |    2     2
--R   (27)  - \|- q2  - q1  + 1 %i i + q2 %i j + q1 %i k
--R                               Type: Quaternion(Complex(Expression(Integer)))
--E 27

———-

5 Assumptions about Variables

Now, in the equations that will result, we will obtain subexpressions of the form \sqrt{θ^{2}}. The positive root will always be appropriate in this formula and one way to get it (in other computer algebra systems) would be to restrict θ to positive values so as to avoid branches that would otherwise have to be properly selected to get an equality. Thus, one means to get the desired simplification is to perform the restriction of θ to positive values. My current understanding is: Axiom does not support variables with properties in the sense that Maple does (the *assume* facility).

6 Problems

Returning to the case at hand, let’s expand a test of equality to the four biquat components and seek a way to prove that each component is zero.

— * —

--S 28 of 43
theta:=_\theta

--R
--R
--R   (28)  \theta
--R                                                       Type: Variable(\theta)
--E 28
--S 29 of 43
testqeq:=map(simplify,siglist(rot(theta,qnv)-exp((-theta/2)*%i*qnv)))_
         ::List Expression Complex Integer
--R
--R
--R   (29)
--R           +-------+
--R           |      2
--R          \|\theta          \theta
--R   [- cos(----------) + cos(------),
--R               2               2
--R                        +-------+
--R          +-------+     |      2
--R          |      2     \|\theta                       \theta
--R    %i q1\|\theta  sin(----------) - %i \theta q1 sin(------)
--R                            2                            2
--R    ---------------------------------------------------------,
--R                              \theta
--R                        +-------+
--R          +-------+     |      2
--R          |      2     \|\theta                       \theta
--R    %i q2\|\theta  sin(----------) - %i \theta q2 sin(------)
--R                            2                            2
--R    ---------------------------------------------------------,
--R                              \theta
--R
--R                                            +-------+
--R            +---------------+ +-------+     |      2
--R            |    2     2      |      2     \|\theta
--R         %i\|- q2  - q1  + 1 \|\theta  sin(----------)
--R                                                2
--R       +
--R                                 +---------------+
--R                         \theta  |    2     2
--R         - %i \theta sin(------)\|- q2  - q1  + 1
--R                            2
--R    /
--R       \theta
--R     ]
--R                                     Type: List(Expression(Complex(Integer)))
--E 29

———- As is quite obvious from this result, the two equations would have been shown to be equal if positivity of the theta variable was enforced. This raises two questions

1. How would one perform a proof such as the above within the allowed operations of Axiom?

2. How can Axiom be called a symbolic algebra system if basic notions of algebraic proof such as restrictions of a variable to the semiring of non-negative integers is disallowed?

How are symbolic manipulations, proofs and solutions to be carried out without this basic notion? I am let to think that Axiom can only be considered a semi-numerical software system.

7 Rules and Pattern Matching

Let’s try another way to answer question 1 by using symbolic pattern matching “rules”. Define a rule that recognizes the specific case of the variable theta appearing in a construct that should be simplified given its positivity.

— * —

--S 30 of 43
posthetaRule:=rule sqrt(theta^2)==theta
--R
--R
--R          +------+
--R          |     2
--R   (30)  \|theta   == theta
--R                       Type: RewriteRule(Integer,Integer,Expression(Integer))
--E 30

———-

But, now we encounter two new problems. First, apparently “map” does not work with rules, so I can’t apply this rule in a simple fashion to the entire List generated by our test case above for any x.

The “map” function only maps functions but we can easily turn a rule into an anonymous function.

— * —

--S 31 of 43
map(x+->posthetaRule(x), [0,sqrt(theta^2),0,sqrt(theta^2)])
--R
--R
--R   (31)  [0,\theta,0,\theta]
--R                                              Type: List(Expression(Integer))
--E 31

———-

But even if we were to split our list above, and try to apply this rule to a component, we get this disappointing result:

— * —

--S 32 of 43
posthetaRule testqeq.1
--R
--R   There are no library operations named posthetaRule
--R      Use HyperDoc Browse or issue
--R                            )what op posthetaRule
--R      to learn if there is any operation containing " posthetaRule " in
--R      its name.
--R
--RDaly Bug
--R   Cannot find a definition or applicable library operation named
--R      posthetaRule with argument type(s)
--R                        Expression(Complex(Integer))
--R
--R      Perhaps you should use "@" to indicate the required return type,
--R      or "$" to specify which version of the function you need.
--E 32

———-

Apparently rules have builtin assumptions about the construction of the expression. If I override this for a real element of the list we are testing I can successfully get one answer.

— * —

--S 33 of 43
[posthetaRule (testqeq.i::Expression Integer) for i in 1..1]
--R
--R
--R   (32)  [0]
--R                                              Type: List(Expression(Integer))
--E 33

———-

But this doesn’t extend to the remainder of the list.

— * —

--S 34 of 43
[posthetaRule (testqeq.i::Expression Integer) for i in 1..4]
--R
--R
--RDaly Bug
--R   Cannot convert from type Expression(Complex(Integer)) to Expression(

--R      Integer) for value
--R                       +-------+
--R         +-------+     |      2
--R         |      2     \|\theta                       \theta
--R   %i q1\|\theta  sin(----------) - %i \theta q1 sin(------)
--R                           2                            2
--R   ---------------------------------------------------------
--R                             \theta
--R
--E 34

———-

because obviously we can’t coerce an actually complex value into the Expression Integer type that the rule seems to want.

Here is one way to create the necessary rule. Like everything else in Axiom rules are objects in some domain. We need to specify a domain for the rule that is compatible with the domain we wish to manipulate. So let’s find out more about this domain:

— * —

--S 35 of 43
)show RewriteRule
--R
--I RewriteRule(Base: SetCategory,R: Join(Ring,...
--R Abbreviation for RewriteRule is RULE
--R This constructor is exposed in this frame.
--R Issue )edit bookvol10.3.pamphlet to see algebra source code for RULE
--R
--R------------------------------- Operations --------------------------------
--R ?=? : (%,%) -> Boolean                coerce : Equation(F) -> %
--R coerce : % -> OutputForm              elt : (%,F,PositiveInteger) -> F
--R ?.? : (%,F) -> F                      hash : % -> SingleInteger
--R latex : % -> String                   lhs : % -> F
--R pattern : % -> Pattern(Base)          retract : % -> Equation(F)
--R rhs : % -> F                          rule : (F,F,List(Symbol)) -> %
--R rule : (F,F) -> %                     ?~=? : (%,%) -> Boolean
--R quotedOperators : % -> List(Symbol)
--R retractIfCan : % -> Union(Equation(F),"failed")
--R suchThat : (%,List(Symbol),(List(F) -> Boolean)) -> %
--R
--E 35

———-

We can provide satisfactory parameters of this domain constructor as follows: — * —

--S 36 of 43

Complex Integer has PatternMatchable Integer
--R
--R
--R   (33)  true
--R                                                                Type: Boolean
--E 36
--S 37 of 43
Expression Complex Integer has FunctionSpace Complex Integer
--R
--R
--R   (34)  true
--R                                                                Type: Boolean
--E 37
--S 38 of 43
Expression Complex Integer has PatternMatchable Integer
--R
--R
--R   (35)  true
--R                                                                Type: Boolean
--E 38

———-

So the following simplification works as hoped:

— * —

--S 39 of 43
posxRule:=(rule sqrt('theta^2)==theta)$RewriteRule(Integer,Complex Integer,_
            Expression Complex Integer)
--R
--R
--R          +------+
--R          |     2
--R   (36)  \|theta   == theta
--R     Type: RewriteRule(Integer,Complex(Integer),Expression(Complex(Integer)))
--E 39
--S 40 of 43
map(x+->posxRule x, testqeq)
--R
--R
--R   (37)  [0,0,0,0]
--R                                     Type: List(Expression(Complex(Integer)))
--E 40

———-

8 Algebraic Domain of Computation

Recall that Axiom does things in a fundamentally algebraic rather than symbolic way. So although the pattern matching approach works, it goes somewhat against the overall philosophy of Axiom. It is possible to do almost anything by using rewrite rules including things that are mathematically incorrect. In this case the rule is simple and the results obvious but if it gets more complicated this can be a problem.

Is there some “axiomatic”, i.e. algebraic way to express the fact that θ is a positive real? well, what we really need is that

\sqrt{θ^{2}} = θ

be satisfied. Now consider that by definition √:

— * —

--S 41 of 43
test (sqrt(x)^2=x)
--R
--R
--R   (38)  true
--R                                                                Type: Boolean
--E 41

———-

and so — * —

--S 42 of 43
test (sqrt(sqrt(x)^2)=sqrt(x))
--R
--R
--R   (39)  true
--R                                                                Type: Boolean
--E 42

———-

for any x. So the expression \sqrt{(x)} is one such thing that has the desired property in the Expression domain. With exactly the appropriate loss of generality, let us just suppose that θ is \sqrt{(β)} for some β. Then if necessary we can also perform a change of variable back to θ.

— * —

--S 43 of 43
eval(eval(testqeq,theta=sqrt(beta)),sqrt(beta)=theta)
--R
--R

--R   (40)  [0,0,0,0]
--R                                     Type: List(Expression(Complex(Integer)))
--E 43
)spool
)lisp (bye)

———-

References^{[1] nothing}

19