Physics - New Features in Maple 2019 - Maplesoft

What's New in Maple 2019

Physics



Maple provides a state-of-the-art environment for algebraic computations in Physics, with emphasis on ensuring that the computational experience is as natural as possible. The theme of the Physics project for Maple 2019 has been the consolidation of the functionality introduced in previous releases, together with significant enhancements to further strengthen the functionality mainly in three areas: 

  • Quantum Mechanics: coherent states, tensor products of states, taylor series of expressions involving anticommutative variables and functions, and several improvements in the normalization and simplification of Commutator and AntiCommutator algebra rules.
  • Tensor computations in general, making Maple 2019 unmatched in the field, covering classical and quantum mechanics, and special and general relativity, using natural tensor input notation and textbook-like display of results. The functionality for tensors is tightly integrated with the full Maple computation system and extensively documented in A Complete Guide for performing tensor computations using Physics.
  • Documentation: besides the new guide for tensor computations, two other new pages, linked in all the help pages of Physics commands, are
  • Physics Updates organizes and presents in one place all formerly scattered links to updates and presentations with examples on the use of Physics.
  • Mini-Course: Computer Algebra for Physicists, is a course that can be used as a tutorial, with 10 sections to be covered in 5 hands-on guided experiences of 2 hours each. The first part, 5 sections, is about Maple 101, while the remaining 5 sections is all about using the Physics package.
 

Overall, the enhancements throughout the entire package increase robustness, versatility and functionality, extending furthermore the range of Physics-related algebraic computations that can be done naturally in a worksheet. The presentation below illustrates both the novelties and the kind of mathematical formulations that can now be performed.
 

As part of its commitment to providing the best possible environment for algebraic computations in Physics, Maplesoft launched a Maple Physics: Research and Development website with Maple 18, which enabled users to download research versions, ask questions, and provide feedback. The results from this accelerated exchange with people around the world have been incorporated into the Physics package in Maple 2019.  

 

Tensor product of Quantum States using Dirac's Bra-Ket Notation 

 Tensor products of Hilbert spaces and related quantum states are relevant in a myriad of situations in quantum mechanics, and in particular regarding quantum information. Tensor products are key in the mathematical formulation of entanglement. Below is a presentation of the design and implementation introduced in Maple 2019, with input/output and examples, organized in four sections: 

  • The basic ideas and design
  • Tensor product notation and the hideketlabel option
  • Entangled States and the Bell basis
  • Entangled States, Operators and Projectors
 

 

References 

[1] Cohen-Tannoudji, C.; Diu, B.; and Laloe,F. Quantum Mechanics, Chapter 2, section F. 

[2] Griffiths, Robert B. Hilbert Space Quantum Mechanics. Quantum Computation and Quantum Information Theory Course, Physics Department, Carnegie Mellon University, 2014. 

 

The basic ideas and design implemented in Maple 2019 

 

Suppose A and B are quantum operators and Ket(A, n), et(B, m) are, respectively, their eigenkets. The following works since the introduction of the Physics package in Maple 

> with(Physics); -1
 

> Setup(op = {A, B})
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[quantumoperators = {A, B}] (1)
 

> `*`(A, `*`(Ket(A, alpha))) = Typesetting:-delayDotProduct(A, Ket(A, alpha))
 

Typesetting:-mprintslash([Physics:-`*`(A, Physics:-Ket(A, alpha)) = `*`(alpha, `*`(Physics:-Ket(A, alpha)))], [Physics:-`*`(A, Physics:-Ket(A, alpha)) = `*`(alpha, `*`(Physics:-Ket(A, alpha)))]) (2)
 

> `*`(B, `*`(Ket(B, beta))) = Typesetting:-delayDotProduct(B, Ket(B, beta))
 

Typesetting:-mprintslash([Physics:-`*`(B, Physics:-Ket(B, beta)) = `*`(beta, `*`(Physics:-Ket(B, beta)))], [Physics:-`*`(B, Physics:-Ket(B, beta)) = `*`(beta, `*`(Physics:-Ket(B, beta)))]) (3)
 

In previous Maple releases, all quantum operators are supposed to act on the same Hilbert space. New in Maple 2019: suppose that A and B act on different, disjointed, Hilbert spaces. 

 

1) To represent that situation, a new keyword in Setup, hilbertspaces, is introduced. With it you can indicate the quantum operators that act on a Hilbert space, say as in hilbertspaces = {{A}, {B}} with the meaning that the operator A acts on one Hilbert space while B acts on another one.  

 

The Hilbert space thus has no particular name (as in 1, 2, 3 ...) and is instead identified by the operators that act on it. There can be one or more, and operators acting on one space can act on other spaces too. The disjointedspaces keyword is a synonym for hilbertspaces and hereafter all Hilbert spaces are assumed to be disjointed. 

 

NOTE: noncommutative quantum operators acting on disjointed spaces commute between themselves, so after setting - for instance - `*`(hilbertspaces, `*`(spaces)) = {{A}, {B}}, automatically, A, B become quantum operators satisfying (see comment (ii) on page 156 of ref.[1]) 

 

 

 

2) Product of Kets and Bras that belong to different Hilbert spaces are understood as tensor products satisfying (see footnote on page 154 of ref. [1]): 

 

`⊗`(Ket(A, alpha), Ket(B, beta)) = `⊗`(Ket(B, beta), Ket(A, alpha))  

 

`⊗`(Bra(A, alpha), Ket(B, beta)) = `⊗`(Ket(B, beta), Bra(A, alpha))  

 

while 

`<>`(`*`(Bra(A, alpha), `*`(Ket(A, alpha))), `*`(Bra(A, alpha), `*`(Ket(A, alpha)))) 

 

3) All the operators of one Hilbert space act transparently over operators, Bras, and Kets of other Hilbert spaces. For example 

 

`*`(A, `*`(Ket(B, n))) = `*`(A, `*`(Ket(B, n))) 

and the same for the Dagger of this equation, that is  

`*`(Bra(B, n), `*`(Dagger(A))) = `*`(Bra(B, n), `*`(Dagger(A))) 

 

Hence, when we write the left-hand sides of the two equations above and press enter, they are automatically rewritten (returned) as the right-hand sides. 

 

4) Every other quantum operator, set as such using Setup, and not indicated as acting on any particular Hilbert space, is assumed to act on all spaces. 

 

5) Notation:  

 

  • Tensor products formed with operators, or Bras and Kets belonging to different Hilbert spaces (set as such using Setup and the keyword hilbertspaces), are now displayed with the symbol ⊗ in between, as in `*`(Ket(A, n), `*`(Ket(B, n))) instead of `*`(Ket(A, n), `*`(Ket(B, n))), and `⊗`(A, B) instead of `*`(A, `*`(B)). The product of an operator A of one space and a Ket of another space Ket(B, n) however, is displayed AA, without ⊗.
 

  • A new Setup option hideketlabel , makes all the labels in Kets and Bras to be hidden when displaying Kets, Bras, and Bracket, so when you set it entering Setup(hideketlabel = true),
       
 

  

is displayed as  

Ket(A, m, n, l) 

 

  • This is the notation frequently used when working with angular momentum or in quantum information, where tensor products of Hilbert spaces are used.
 

Design details 

 

The commutativity of the eigenkets of A and B is consistent with , see footnote on page 154 of ref. [1]. 

 

Taking advantage of this commutativity of Bras and Kets belonging to disjointed spaces, during the computer algebra session (on the worksheet) the ordering of their products in the output happens automatically and systematically, it is always the same, and according to the following rules: suppose there are two Hilbert subspaces, then: 

 

  • Within a product, contiguous Kets are grouped per Hilbert subspace
 

  • The Hilbert subspaces are ordered by sorting, alphabetically, the operators that act on that subspace.
 

 

Example: if one Hilbert subspace has operators {A, E} acting on it and the other has operators {B, C} then a product of contiguous eigenkets of these operators is sorted as 

 

`⊗`(`⊗`(`⊗`(Ket(A), Ket(E)), Ket(B)), Ket(C)) 

 

Where the first pair of Kets belong to the first Hilbert subspace and the other pair to the second subspace, and where the first subspace is the one whose operands are alphabetically sorted first than the operands of the second subspace (in this example A is sorted before B) and within a subspace, the Kets are also sorted alphabetically (so A before E, then in the second subspace B before C). 

 

Regarding the notation for the Dagger of a tensor product of states, say Dagger(`*`(Ket(A, alpha), Ket(B, beta))) the standard convention for tensor products is to preserve the order, as in Dagger(`*`(Ket(A, alpha), Ket(B, beta))) = `⊗`(Bra(A, alpha), Bra(B, beta))representing an exception to the “reverse the order” rule of the Dagger operation. This is conventional, in that Kets and Bras belonging to disjointed spaces actually commute. This convention, however is notationally important for two reasons 

 

  • 1. Frequently, tensor product of states are written not as a product of Kets but as a single ket with many quantum numbers (multidimensional Hilbert space), for example as Ket(C, alpha, beta), and Dagger(Ket(C, alpha, beta))as Bra(C, alpha, beta), that is, preserving the order of the quantum numbers, the first index refers to the first Hilbert subspace and the second index to the other one, in both the Ket and the Bra (its Dagger). So, when these multi-index Kets can be expressed as tensor products, the ordering of the Hilbert subspaces is preserved.
 

  • 2. This preservation of the ordering of the Hilbert subspaces in a tensor product of Kets is also relevant in connection with the new Setup keyword hideketlabel, which requires the ordering of subspaces to be preserved in order to have non-ambiguous notation when the quantum numbers of Kets belonging to identical subspaces (e.g. qubits) are the same.
 

 

Example: we know that 

 

 

In the left-hand and right-hand sides of the expression above, the ordering of the Hilbert subspaces is not the same. If we now omit the labels A and B, we would have  

 

`<>`(`⊗`(Ket(0), Ket(1)), `⊗`(et(0), Ket(1))) 

which would be misleading. Likewise,  

`⊗`(Ket(A, 0), Ket(B, 1)) = `⊗`(Ket(B, 1), Ket(A, 0)) 

and removing the labels we would get the misleading  

`⊗`(Ket(0), Ket(1)) = `⊗`(Ket(0), Ket(1)) 

 

From all this we see that, in order to make sense of the notation without labels, it is necessary to preserve the ordering of the Hilbert subspaces present in a tensor product, and also when taking the Dagger of a Ket. Accordingly, within tensor products, for instance in these examples, the system will always write Kets of the subspace A before Kets of the subspace B. 

 

Tensor product notation and the hideketlabel option 

 

According to the design section, set now two disjointed Hilbert spaces with operators A, C acting on one of them and B, C on the other one (you can think of  C = `⊗`(A, B)) 

 

> Setup(hilbertspaces = {{A, C}, {B, C}})
 

[disjointedspaces = {{A, C}, {B, C}}] (4)
 

 

Consider a tensor product of Kets, each of which belongs to one of these different spaces, note the new notation using 

> `*`(Ket(A, 1), `*`(Ket(B, 0)))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0))], [Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0))]) (5)
 

  • As explained in the Details of the design section, the ordering of the Hilbert spaces in tensor products is now preserved: Bras (Kets) of the first space always appear before Bras (Kets) of the second space. For example, construct a projector into the state `*`(Ket(A, 1), `*`(Ket(B, 0)))
 

> `*`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)), `*`(Dagger(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)))))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))], [Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physi... (6)
 

You see that in the product of Bras, and also in the product of Kets, A comes first, then B. 


Remark
: some textbooks prefer a dyadic style for sorting the operands in products of Bras and Kets that belong to different spaces, for example, `⊗`(`*`(Ket(A, 1), `*`(Bra(A, 1))), `⊗`(Ket(B, 0), Bra(B, 0))) instead of the projector sorting style of `*`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)), `*`(Dagger(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0))))). Both reorderings of Kets and Bras are mathematically equal. 

 

  • Because that ordering is preserved, one can now hide the label of Bras and Kets without ambiguity, as it is usual in textbooks (e.g. in Quantum Information). For that purpose use the new keyword option hideketlabel
 

> Setup(hide = true)
 

 

 

`*`(`* Partial match of  '`, `*`(hide, `*`(`' against keyword '`, `*`(hideketlabel, `*`(`' `)))))
_______________________________________________________
[hideketlabel = true] (7)
 

The display for `*`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)), `*`(Dagger(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0))))) is now 

> Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))], [Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physi... (8)
 

Important: this new option only hides the label while displaying the Bra or Ket. The label, however, is still there, both in the input and in the output. One can "see" what is behind this new display using show, that works the same way as it does in the context of CompactDisplay. The actual contents being displayed in Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0)) is thus `*`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)), `*`(Dagger(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0))))) 

> show
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))], [Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physi... (9)
 

Operators of each of these spaces act on their eigenkets as usual. Here we distribute over both sides of an equation, using `*` on the left-hand side, to see the product uncomputed, and `.` on the right-hand side to see it computed: 

> (`*` = `.`)(A, Ket(A, 1)); 1
 

Typesetting:-mprintslash([Physics:-`*`(A, Physics:-Ket(A, 1)) = Physics:-Ket(A, 1)], [Physics:-`*`(A, Physics:-Ket(A, 1)) = Physics:-Ket(A, 1)]) (10)
 

> (`*` = `.`)(A, Ket(A, 0)); 1
 

Typesetting:-mprintslash([Physics:-`*`(A, Physics:-Ket(A, 0)) = 0], [Physics:-`*`(A, Physics:-Ket(A, 0)) = 0]) (11)
 

  • The tensor product of operators belonging to different Hilbert spaces is also displayed using ⊗
 

> `*`(A, `*`(B))
 

Typesetting:-mprintslash([Physics:-`*`(A, B)], [Physics:-`*`(A, B)]) (12)
 

  • As mentioned in the preceding design section, using the commutativity between operators, Bras and Kets that belong to different Hilbert spaces, within a product, operators are placed contiguous to the Kets and Bras belonging to the space where the operator acts. For example, consider the delayed product represented using the start `*` operator
 

> '`*`(Physics:-`*`(A, B), `*`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))))'
 

Typesetting:-mprintslash([Physics:-`*`(A, B, Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))], [Physics:-`*`(A, B, Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(... (13)
 

Release the product 

> %
 

Typesetting:-mprintslash([Physics:-`*`(A, Physics:-Ket(A, 1), B, Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0))], [Physics:-`*`(A, Physics:-Ket(A, 1), B, Physics:-Ket(B, 0), Physics:-Bra(... (14)
 

The same operation but now using the dot product `.` operator. Start by delaying the operation 

> 'Typesetting:-delayDotProduct(Physics:-`*`(A, B), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0)))'
 

Typesetting:-mprintslash([`.`(Physics:-`*`(A, B), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0)))], [Typesetting:-delayDotProduct(Physics:-`*`(A, B), Phys... (15)
 

Recalling that this product is mathematically the same as %, and that  

> Typesetting:-delayDotProduct(B, Ket(B, 0))
 

0 (16)
 

by releasing the delayed product 'Typesetting:-delayDotProduct(Physics:-`*`(A, B), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0)))' we have 

> Typesetting:-delayDotProduct(Physics:-`*`(A, B), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0), Physics:-Bra(A, 1), Physics:-Bra(B, 0)))
 

0 (17)
 

Reset hideketlabel 

> Setup(hideketlabel = false)
 

[hideketlabel = false] (18)
 

Implementation details 

 

  • When defining a disjointed Hilbert space that contains operators belonging to one of the previously defined disjointed spaces, if the previously defined space is a subset of the one being defined, the former definition is removed and only the latter is kept.
 

  • There are three new routines related to operators acting on disjointed spaces
 

> Library:-BelongToDifferentDisjointedSpaces(A, B)
 

true (19)
 

  • and, in connection with the possibility of indicating that two operators act on the same disjointed space, there is a new routine to tell when the operators belong to (act on) the same space
 

> Library:-BelongToTheSameDisjointedSpace(A, B)
 

false (20)
 

> Library:-BelongToTheSameDisjointedSpace(A, C)
 

true (21)
 

  • Finally, there is a routine to tell whether a sequence of objects belong (all of them) to (the same or different) disjointed spaces, i.e. where set as such using Setup and the hilbertpaces keyword
 

> Library:-BelongToDisjointedSpaces(A, B, C)
 

true (22)
 

> Setup(op = E)
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[quantumoperators = {A, B, C, E}] (23)
 

> Library:-BelongToDisjointedSpaces(E)
 

false (24)
 

> Library:-BelongToDisjointedSpaces(A, B, E)
 

false (25)
 

These three Physics:-Library routines, are the ones used internally by the Physics package to make decisions.
 

  • Coordinates of a coordinate system defined using the Coordinates, or equivalently the Setup, command, when they are also set to as quantum operators, are automatically assumed to act on disjointed spaces. For example,
 

> Setup(dimension = 3, signature = `+`, coordinates = cartesian, quantumoperators = {X}, quiet)
 

[coordinatesystems = {X}, dimension = 3, quantumoperators = {A, B, C, E, x, y, z}, signature = `+ + +`] (26)
 

> Library:-BelongToDifferentDisjointedSpaces(x, y)
 

true (27)
 

  • Remove the coordinates from the quantum operators set
 

> Setup(clear, op = {X})
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[quantumoperators = {A, B, C, E}] (28)
 

  • The following swapped products on the left-hand and right-hand sides (with their evaluation delayed due to enclosing the equation between ' ' ) are automatically reordered in one and the same way when we remove the delay, as explained in the details of the design section.
 

 

Example 

 

> Setup(hide = false)
 

 

 

`*`(`* Partial match of  '`, `*`(hide, `*`(`' against keyword '`, `*`(hideketlabel, `*`(`' `)))))
_______________________________________________________
[hideketlabel = false] (29)
 

> '`*`(Ket(A, n), `*`(Ket(B, m))) = `*`(Ket(A, n), `*`(Ket(B, m)))'
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:-`*`(Physics:-Ket(B, m), Physics:-Ket(A, n))], [Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:... (30)
 

Remove now the evaluation delay and the ordering on the right-hand side is automatically rearranged as in the left-hand side. 

> %
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m))], [Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:... (31)
 

The same using the dot operator `.` 

> 'Typesetting:-delayDotProduct(Ket(A, n), Ket(B, m)) = Typesetting:-delayDotProduct(Ket(B, m), Ket(A, n))'
 

Typesetting:-mprintslash([`.`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = `.`(Physics:-Ket(B, m), Physics:-Ket(A, n))], [Typesetting:-delayDotProduct(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Typesettin... (32)
 

> Typesetting:-delayDotProduct(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Typesetting:-delayDotProduct(Physics:-Ket(B, m), Physics:-Ket(A, n))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m))], [Physics:-`*`(Physics:-Ket(A, n), Physics:-Ket(B, m)) = Physics:... (33)
 

NOTE: the dot product operator, `.`, is used to perform contractions or attachments in the space of quantum states. Therefore, in the case of tensor products, it returns using the star product operator `*`, in that there is no meaning for the contraction of tensors of different (disjointed) Hilbert spaces. 

 

Regarding the product of a Bra and a Ket belonging to disjointed spaces, we also have, automatically,  

> '`*`(Bra(A, n), `*`(Ket(B, n))) = `*`(Bra(A, n), `*`(Ket(B, n)))'
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(A, n), Physics:-Ket(B, n)) = Physics:-`*`(Physics:-Ket(B, n), Physics:-Bra(A, n))], [Physics:-`*`(Physics:-Bra(A, n), Physics:-Ket(B, n)) = Physics:... (34)
 

> %
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(B, n), Physics:-Bra(A, n)) = Physics:-`*`(Physics:-Ket(B, n), Physics:-Bra(A, n))], [Physics:-`*`(Physics:-Ket(B, n), Physics:-Bra(A, n)) = Physics:... (35)
 

So the left-hand side is rewritten as the right-hand side, and `⊗`(Bra(A, n), Ket(B, n))is not a "scalar product", but an operator in the tensor product of spaces, since A and B belong to different disjointed spaces. 

 

Enclose again the input with ' ' to delay its evaluation 

> '`*`(A, `*`(Ket(B, n))) = `*`(A, `*`(Ket(B, n)))'
 

Typesetting:-mprintslash([Physics:-`*`(A, Physics:-Ket(B, n)) = Physics:-`*`(Physics:-Ket(B, n), A)], [Physics:-`*`(A, Physics:-Ket(B, n)) = Physics:-`*`(Physics:-Ket(B, n), A)]) (36)
 

Release the evaluation 

> %
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(B, n), A) = Physics:-`*`(Physics:-Ket(B, n), A)], [Physics:-`*`(Physics:-Ket(B, n), A) = Physics:-`*`(Physics:-Ket(B, n), A)]) (37)
 

In the output above we see that Typesetting:-delayDotProduct(A, Ket(B, n)) is not interpreted as contraction between an operator and Ket, but as the product of A = `*`(A, `*`(`𝕀`[B])) acting on Ket(B, n) where `𝕀`[B] is the identity (projector) onto the B space. That is, an operator of one disjointed space acts transparently over a Bra or a Ket of a different disjointed space. The same happens with Dagger(A) just that, while A moves to the right, jumping over a Bra or Ket (see %), Dagger(A) moves to the left: 

> '`*`(Dagger(A), `*`(Bra(B, n))) = `*`(Dagger(A), `*`(Bra(B, n)))'
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (38)
 

> %
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (39)
 

NOTE: Although determining "who is the Dagger of who" is arbitrary, this implementation follows what we do with paper and pencil: operators act to their right while those having an explicit Dagger act to their left.  

 

Finally, the notation used for tensor products of operators is the same one used for tensor products of Bras and Kets: 

> `*`(A, `*`(B))
 

Typesetting:-mprintslash([Physics:-`*`(A, B)], [Physics:-`*`(A, B)]) (40)
 

As explained in the Details of the design section, the ordering of the Hilbert spaces in tensor products is now preserved, so taking the Dagger does not swap the operands in this product: 

> Dagger(Physics:-`*`(A, B))
 

Typesetting:-mrow(Typesetting:-mi( (41)
 

Entangled States and the Bell basis 

 

With the introduction of disjointed Hilbert spaces in Maple 2019 it is possible to represent entangled quantum states in a simple way, basically as done with paper and pencil.  

 

Recalling the Hilbert spaces set at this point are, 

> Setup(hilbert)
 

 

 

`*`(`* Partial match of  '`, `*`(hilbert, `*`(`' against keyword '`, `*`(hilbertspaces, `*`(`' `)))))
_______________________________________________________
[disjointedspaces = {{A, C}, {B, C}}] (42)
 

where C acts on the tensor product of the spaces where A and B act. A state of C can then always be written as 

> Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Ket(A, j), `*`(Ket(B, p)))), j), p)
 

Typesetting:-mprintslash([Physics:-Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j), p)], [Physics:-Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Physics:-`... (43)
 

where M[j, p] is a matrix of complex coefficients. Bra states of C are formed as usual taking the Dagger 

> Dagger(Physics:-Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j), p))
 

Typesetting:-mprintslash([Physics:-Bra(C, m, n) = Sum(Sum(`*`(conjugate(M[j, p]), `*`(Physics:-`*`(Physics:-Bra(A, j), Physics:-Bra(B, p)))), j), p)], [Physics:-Bra(C, m, n) = Sum(Sum(`*`(conjugate(M[... (44)
 

 

  • By definition, all states Ket(C, alpha, beta) that can be written exactly as `⊗`(Ket(A, alpha), Ket(B, beta)), that is, the product of an arbitrary state of the subspace A and another of the subspace B, are product states, and all the other ones are entangled states. Entanglement is a property that is independent of the basis `⊗`(Ket(A, j), Ket(B, p))used in Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Ket(A, j), `*`(Ket(B, p)))), j), p).

    The physical interpretation is the standard one: when the state of a system constituted by two subsystems A and B is represented by a
    product state, the properties of the subsystem A are well defined and all given by while those for the subsystem B by . When the system is in an entangled state one typically cannot assign definite properties to the individual subsystems A or B, each subsystem has no independent reality.

    To determine whether a state Ket(C, alpha, beta) is entangled it then suffices to check the rank R of the matrix M[j, p] : when R = 1 the state is a
    product state, otherwise it is an entangled state. When the state being analyzed belongs to the tensor product of two subspaces, R = 1 is equivalent to having the determinant of M[j, p] equal to 0. The condition R = 1, however, is more general, and suffices to determine whether a state is a product state also on a Hilbert space that is the tensor product of three or more subspaces: `ℋ` = `⊗`(`⊗`(diff(`ℋ`(x), x), diff(`ℋ`(x), x, x)), diff(`ℋ`(x), x, x, x)) .. diff(`ℋ`(x), [`$`(x, n)]), in which case the matrix M will have more rows and columns and a determinant equal to 0 would only warrant the possibility of factorizing one Ket.
 

 

Example: the Bell basis for a system of two qubits 

 

Consider a 2-dimensional space of states acted upon by the operator A, and let B act upon another, disjointed, Hilbert space that is a replica of the Hilbert space on which A acts. Set the dimensions of A, B and C respectively equal to 2, 2 and 2x2  

> Setup(quantumbasisdimension = {A = 2, B = 2, C[1] = 2, C[2] = 2})
 

[quantumbasisdimension = {A = 2, B = 2, C[1] = 2, C[2] = 2}] (45)
 

The system C with the two subsystems A and B represents a two qubits system. The standard basis for C can be constructed in a natural way from the basis of Kets of A and B, {Ket(A, 0), Ket(A, 1), Ket(B, 0), Ket(B, 1)}, by taking their tensor products: 

> seq(seq(`*`(Ket(A, j), `*`(Ket(B, k))), k = 0 .. 1), j = 0 .. 1)
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)), Physics:-`*`... (46)
 

Set a more mathematical display for the imaginary unit 

> interface(imaginaryunit = i); -1
 

 

The four entangled Bell states also form a basis of C and are given by 

> Setup(op = `ℬ`)
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[quantumoperators = {`ℬ`, A, B, C, E}] (47)
 

> Ket(`ℬ`, 0) = `/`(`*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 0))), `*`(Ket(A, 1), `*`(Ket(B, 1))))), `*`(('sqrt')(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 0) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1)))), `*`(sqrt(2)))], [Physics... (48)
 

> Ket(`ℬ`, 1) = `/`(`*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 1))), `*`(Ket(A, 1), `*`(Ket(B, 0))))), `*`(('sqrt')(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 1) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)))), `*`(sqrt(2)))], [Physics... (49)
 

> Ket(`ℬ`, 2) = `/`(`*`(i, `*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 1))), `-`(`*`(Ket(A, 1), `*`(Ket(B, 0))))))), `*`(('sqrt')(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 2) = `/`(`*`(I, `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), `-`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)))))), `*`(sqrt(2)... (50)
 

> Ket(`ℬ`, 3) = `/`(`*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 0))), `-`(`*`(Ket(A, 1), `*`(Ket(B, 1)))))), `*`(('sqrt')(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 3) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), `-`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1))))), `*`(sqrt(2)))], [Ph... (51)
 

There is no standard convention for the four linear combinations of the right-hand sides above defining the Bell states. The convention used here relates to the definition of these states using the Pauli matrices as shown further below. Regardless of the convention used, the Bell basis is orthonormal. That can be verified by taking dot products, for example: 

> Typesetting:-delayDotProduct(Dagger(Physics:-Ket(`ℬ`, 0) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1)))), `*`(sqrt(2))))...
 

1 = 1 (52)
 

In steps, perform the same operation but using the star (`*`) operator, so that the contraction is represented but not performed 

> `*`(Dagger(Physics:-Ket(`ℬ`, 0) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1)))))))...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`ℬ`, 0), Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(Physics:-`*`(`+`(Physics:-`*`(Physics:-Bra(A, 0), Physics:-Bra(B, 0)), Physics:-`*... (53)
 

Evaluate now the result at `*` = `.`, that is transforming the star product into a dot product 

> eval(Physics:-`*`(Physics:-Bra(`ℬ`, 0), Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(Physics:-`*`(`+`(Physics:-`*`(Physics:-Bra(A, 0), Physics:-Bra(B, 0)), Physics:-`*`(Physics:-Bra(A, 1),...
 

1 = 1 (54)
 

> `*`(Dagger(Physics:-Ket(`ℬ`, 0) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1)))))))...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`ℬ`, 0), Physics:-Ket(`ℬ`, 1)) = `+`(`*`(`/`(1, 2), `*`(Physics:-`*`(`+`(Physics:-`*`(Physics:-Bra(A, 0), Physics:-Bra(B, 0)), Physics:-`*... (55)
 

> eval(Physics:-`*`(Physics:-Bra(`ℬ`, 0), Physics:-Ket(`ℬ`, 1)) = `+`(`*`(`/`(1, 2), `*`(Physics:-`*`(`+`(Physics:-`*`(Physics:-Bra(A, 0), Physics:-Bra(B, 0)), Physics:-`*`(Physics:-Bra(A, 1),...
 

0 = 0 (56)
 

The Bell basis and its relation with the Pauli matrices 

 

The Bell basis can be constructed departing from Ket(`ℬ`, 0) using the Pauli matrices sigma[j]. For that purpose, using a Vector representation for Ket(A, j),  

> Ket(B, 0) = Vector([1, 0]), Ket(B, 1) = Vector([0, 1])
 

Typesetting:-mprintslash([Physics:-Ket(B, 0) = Vector[column]([[Typesetting:-mn( (57)
 

Multiplying Ket(B, 0)by each of the sigma[j] Pauli matrices and performing the matrix operations we have 

>
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (58)
 

>
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (59)
 

In this result we see that sigma[1] and sigma[2] flip the state, transforming Ket(B, 0) into Ket(B, 1), sigma[2] also multiplies by the imaginary unit I, while sigma[3] leaves the state Ket(B, 0)unchanged.  

We can express all that by removing from the Vector representations shown in Ket(B, 0) = Vector([1, 0]), Ket(B, 1) = Vector([0, 1]). For that purpose, create a list of substitution equations 

>
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (60)
 

The action of sigma[j] in Ket(B, 0)is then given by 

>
 

Typesetting:-mprintslash([[Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`... (61)
 

For , performing the same steps, the action of the Pauli matrices on it is 

>
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (62)
 

>
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (63)
 

>
 

Typesetting:-mprintslash([[Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 1)) = Physics:-Ket(B, 0), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 1)) = `+`(`-`(`*`(`+`(I), `*`(Physics:-Ket(B, 0))... (64)
 

To obtain the other three Bell states using the results and , indicate to the system that the Pauli matrices operate in the subspace where B operates 

> Setup(hilbert = {{B, C, Psigma}})
 

 

 

`*`(`* Partial match of  '`, `*`(hilbert, `*`(`' against keyword '`, `*`(hilbertspaces, `*`(`' `)))))
_______________________________________________________
[disjointedspaces = {{A, C}, {B, C, Physics:-Psigma}}] (65)
 

 

Multiplying Ket(`ℬ`, 0) given in Ket(`ℬ`, 0) = `/`(`*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 0))), `*`(Ket(A, 1), `*`(Ket(B, 1))))), `*`(('sqrt')(2))) by each of the three sigma[j] we get the other three Bell states 

> Physics:-Ket(`ℬ`, 0) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1)))), `*`(sqrt(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 0) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-... (66)
 

> `*`(Psigma[1], `*`(Physics:-Ket(`ℬ`, 0))) = `+`(`*`(`/`(1, 2), `*`(Psigma[1], `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1),...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(Physics:-`*`(Physics:-Psigma[1], `+`(Physics:-`*`(Physics:-Ket(A, ... (67)
 

Substitute in this result the first equations of and  

> [Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`*`(Physics:-Psigma[3], Phy...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 1)) = Physics:-Ket(B, 0)], [Physics:-`*`(Physics:-... (68)
 

> map(rhs = lhs, [Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 1)) = Physics:-Ket(B, 0)])
 

Typesetting:-mprintslash([[Physics:-Ket(B, 1) = Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)), Physics:-Ket(B, 0) = Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 1))]], [[Physics:-Ket(B, 1) ... (69)
 

> subs([Physics:-Ket(B, 1) = Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)), Physics:-Ket(B, 0) = Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 1))], Physics:-`*`(Physics:-Psigma[1], Physics:-K...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(Physics:-`*`(Physics:-Psigma[1], `+`(Physics:-`*`(Physics:-Ket(A, ... (70)
 

> factor(Simplify(Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(Physics:-`*`(Physics:-Psigma[1], `+`(Physics:-`*`(Physics:-Ket(A, 0), Physic...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), Physics:... (71)
 

This is Ket(`ℬ`, 1) defined in Ket(`ℬ`, 1) = `/`(`*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 1))), `*`(Ket(A, 1), `*`(Ket(B, 0))))), `*`(('sqrt')(2))) 

> Physics:-Ket(`ℬ`, 1) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)))), `*`(sqrt(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 1) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), Physics:-`*`(Physics:-Ket(A, 1), Physics:-... (72)
 

> `+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)), `-`(Physics:-Ket(`ℬ`, 1))) = 0
 

Typesetting:-mprintslash([`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)), `-`(Physics:-Ket(`ℬ`, 1))) = 0], [`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(`ℬ`, 0)), `-`(P... (73)
 

Multiplying now by sigma[2] and substituting Ket(B, j) using the `^`(2, nd) equations of and we get Ket(`ℬ`, 1) 

> `*`(Psigma[2], `*`(Physics:-Ket(`ℬ`, 0))) = `+`(`*`(`/`(1, 2), `*`(Psigma[2], `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1),...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(Physics:-`*`(Physics:-Psigma[2], `+`(Physics:-`*`(Physics:-Ket(A, ... (74)
 

> [Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`*`(Physics:-Psigma[3], Phy...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 1)) = `+`(`-`(`*`(`+`(I), `*`(Physics... (75)
 

> zip(isolate, [Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 1)) = `+`(`-`(`*`(`+`(I), `*`(Physics:-Ket(B, 0))...
 

Typesetting:-mprintslash([[Physics:-Ket(B, 1) = `+`(`-`(`*`(`+`(I), `*`(Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)))))), Physics:-Ket(B, 0) = `*`(I, `*`(Physics:-`*`(Physics:-Psigma[2], Physi... (76)
 

> factor(Simplify(subs([Physics:-Ket(B, 1) = `+`(`-`(`*`(`+`(I), `*`(Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)))))), Physics:-Ket(B, 0) = `*`(I, `*`(Physics:-`*`(Physics:-Psigma[2], Physics:-K...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(`ℬ`, 0)) = `+`(`-`(`*`(`+`(`*`(`/`(1, 2), `*`(I))), `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(Physics:-`*`(Physics:-Ket(A, 0), Phy... (77)
 

The above is Ket(`ℬ`, 2) defined in Ket(`ℬ`, 2) = `/`(`*`(i, `*`(`+`(`*`(Ket(A, 0), `*`(Ket(B, 1))), `-`(`*`(Ket(A, 1), `*`(Ket(B, 0))))))), `*`(('sqrt')(2))) 

> Physics:-Ket(`ℬ`, 2) = `/`(`*`(I, `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), `-`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 0)))))), `*`(sqrt(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 2) = `*`(`*`(`/`(1, 2), `*`(I)), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)), `-`(Physics:-`*`(Physics:-Ket(A, ... (78)
 

> Expand(`+`(Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(`ℬ`, 0)) = `+`(`-`(`*`(`+`(`*`(`/`(1, 2), `*`(I))), `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 1)...
 

Typesetting:-mprintslash([`+`(Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(`ℬ`, 0)), `-`(Physics:-Ket(`ℬ`, 2))) = 0], [`+`(Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(`ℬ`, 0)), `-`(P... (79)
 

Finally, multiplying Ket(`ℬ`, 2) by sigma[3] 

> `*`(Psigma[3], `*`(Physics:-Ket(`ℬ`, 0))) = `+`(`*`(`/`(1, 2), `*`(Psigma[3], `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), Physics:-`*`(Physics:-Ket(A, 1),...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(`ℬ`, 0)) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(Physics:-`*`(Physics:-Psigma[3], `+`(Physics:-`*`(Physics:-Ket(A, ... (80)
 

Substituting  

> [Physics:-`*`(Physics:-Psigma[1], Physics:-Ket(B, 0)) = Physics:-Ket(B, 1), Physics:-`*`(Physics:-Psigma[2], Physics:-Ket(B, 0)) = `*`(I, `*`(Physics:-Ket(B, 1))), Physics:-`*`(Physics:-Psigma[3], Phy...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 0)) = Physics:-Ket(B, 0), Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 1)) = `+`(`-`(Physics:-Ket(B, 1)))], [Physics:-`*`... (81)
 

> (rhs = lhs)((Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 0)) = Physics:-Ket(B, 0), Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 1)) = `+`(`-`(Physics:-Ket(B, 1))))[1]), (rhs = lhs)(`+`(`-`((P...
 

Typesetting:-mprintslash([Physics:-Ket(B, 0) = Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 0)), Physics:-Ket(B, 1) = `+`(`-`(Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 1))))], [Physics:-Ket... (82)
 

We get  

> factor(Simplify(subs(Physics:-Ket(B, 0) = Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 0)), Physics:-Ket(B, 1) = `+`(`-`(Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(B, 1)))), Physics:-`*`(Physic...
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(`ℬ`, 0)) = `+`(`-`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0))),... (83)
 

which is Ket(`ℬ`, 3) 

> Physics:-Ket(`ℬ`, 3) = `/`(`*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), `-`(Physics:-`*`(Physics:-Ket(A, 1), Physics:-Ket(B, 1))))), `*`(sqrt(2)))
 

Typesetting:-mprintslash([Physics:-Ket(`ℬ`, 3) = `+`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0)), `-`(Physics:-`*`(Physics:-Ket(A, 1), Physi... (84)
 

> Expand(`+`(Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(`ℬ`, 0)) = `+`(`-`(`*`(`/`(1, 2), `*`(`^`(2, `/`(1, 2)), `*`(`+`(`-`(Physics:-`*`(Physics:-Ket(A, 0), Physics:-Ket(B, 0))), Physics:-`*`(P...
 

Typesetting:-mprintslash([`+`(Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(`ℬ`, 0)), `-`(Physics:-Ket(`ℬ`, 3))) = 0], [`+`(Physics:-`*`(Physics:-Psigma[3], Physics:-Ket(`ℬ`, 0)), `-`(P... (85)
 

Reset the symbol representing imaginary unit to use i as an index in the next section 

> interface(imaginaryunit = I); -1
 

Entangled States, Operators and Projectors 

 

Consider a fourth operator, H, that is Hermitian and acts on the same space of C, has the same dimension, and `ℋ`, `ℍ` are its mean values in an entangled and product states respectively. 

> Setup(additionally, hermitian = H, basisdimension = {H[1] = 2, H[2] = 2}, hilbertspaces = {{A, C, H}, {B, C, H}}, realobjects = {`ℍ`, `ℋ`})
 

 

 

 

`*`(`* Partial match of  '`, `*`(hermitian, `*`(`' against keyword '`, `*`(hermitianoperators, `*`(`' `)))))
`*`(`* Partial match of  '`, `*`(basisdimension, `*`(`' against keyword '`, `*`(quantumbasisdimension, `*`(`' `)))))
_______________________________________________________
[disjointedspaces = {{A, C, H}, {B, C, H}, {B, C, Physics:-Psigma}}, hermitianoperators = {H}, quantumbasisdimension = {A = 2, B = 2, C[1] = 2, C[2] = 2, H[1] = 2, H[2] = 2}, realobjects = {`ℍ`, ...
[disjointedspaces = {{A, C, H}, {B, C, H}, {B, C, Physics:-Psigma}}, hermitianoperators = {H}, quantumbasisdimension = {A = 2, B = 2, C[1] = 2, C[2] = 2, H[1] = 2, H[2] = 2}, realobjects = {`ℍ`, ...
(86)
 

To operate in a practical way with these operators, Bras and Kets, bracket rules reflecting their relationship are necessary. From the definition of C as acting on the tensor product of spaces where A and B act (see Ket(C, m, n) = Sum(Sum(`*`(M[j, p], `*`(Ket(A, j), `*`(Ket(B, p)))), j), p)) and taking into account the dimensions specified for A, B and C we have 

> Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Ket(A, j), `*`(Ket(B, p)))), j = 0 .. 1), p = 0 .. 1)
 

Typesetting:-mprintslash([Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .. 1), p = 0 .. 1)], [Physics:-Ket(C, a, b) = Sum(Sum(`*`... (87)
 

> Typesetting:-delayDotProduct(Bra(A, k), Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .. 1), p = 0 .. 1))
 

Typesetting:-mprintslash([Physics:-Bracket(Physics:-Bra(A, k), Physics:-Ket(C, a, b)) = Sum(`*`(M[a, k, b, p], `*`(Physics:-Ket(B, p))), p = 0 .. 1)], [Physics:-Bracket(Physics:-Bra(A, k), Physics:-Ke... (88)
 

> Typesetting:-delayDotProduct(Bra(B, k), Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .. 1), p = 0 .. 1))
 

Typesetting:-mprintslash([Physics:-Bracket(Physics:-Bra(B, k), Physics:-Ket(C, a, b)) = Sum(`*`(M[a, j, b, k], `*`(Physics:-Ket(A, j))), j = 0 .. 1)], [Physics:-Bracket(Physics:-Bra(B, k), Physics:-Ke... (89)
 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(A, k), Bra(B, l)), Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .....
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(A, k), Physics:-Bracket(Physics:-Bra(B, l), Physics:-Ket(C, a, b))) = M[a, k, b, l]], [Physics:-`*`(Physics:-Bra(A, k), Physics:-Bracket(Physics:-Br... (90)
 

The bracket rules for A, B and C are the first two of these; Set these rules, so that the system can take them into account 

> Setup(Physics:-Bracket(Physics:-Bra(A, k), Physics:-Ket(C, a, b)) = Sum(`*`(M[a, k, b, p], `*`(Physics:-Ket(B, p))), p = 0 .. 1), Physics:-Bracket(Physics:-Bra(B, k), Physics:-Ket(C, a, b)) = Sum(`*`(...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (91)
 

If we now recompute Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(A, k), Bra(B, l)), Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 ....., the left-hand side is also computed 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(A, k), Bra(B, l)), Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .....
 

M[a, k, b, l] = M[a, k, b, l] (92)
 

Suppose now that you want to compute with the Hermitian operator H, that operates on the same space as C, both using C and the operators A and B, as in 

 

Bracket(Bra(C, I, j), H, Ket(C, k, l)) = `ℋ`[i, j, k, l] 

 

`*`(`⊗`(Bra(A, I), Bra(B, j)), `*`(H, `*`(`⊗`(Ket(A, k), Ket(B, l))))) = `ℍ`[I, j, k, l] 

 

where `ℋ`[i, j, k, l] = `ℍ`[I, j, k, l] when Ket(C, a, b) is a product (not entangled) state.  

 

To compute taking into account Bracket(Bra(C, I, j), H, Ket(C, k, l)) = `ℋ`[I, j, k, l] it suffices to set a bracket rule 

> Setup(%Bracket(Bra(C, a, b), H, Ket(C, c, d)) = `ℋ`[a, b, c, d], real = `ℋ`)
 

 

 

`*`(`* Partial match of  '`, `*`(real, `*`(`' against keyword '`, `*`(realobjects, `*`(`' `)))))
_______________________________________________________
Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(93)
 

After that, 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(C, j, k), H), Ket(C, m, n))
 

`ℋ`[j, k, m, n] (94)
 

Regarding , since H belongs to the tensor product of spaces A and B, it can be an entangled operator, one that you cannot represent just as a product of one operator acting on A times another one acting on B. A computational representation for the operator `*`(Bra(B, j), `*`(H, `*`(Ket(A, k)))) (that is not just itself or as abstract) is not possible in the general case. For that you can use a different feature: define the action of the operator H on Kets of A and B.  

 

Basically, we want: 

 

 

 

A program sketch for that would be: 


if H is applied to a Ket of A or B and it still has not 4 indices then
 

   if H itself is indexed then
       return H with its indices followed by the index of the Ket
   else
 

       return H indexed by the index of the Ket;
otherwise
   return the dot product operation uncomputed, unevaluated
 

 

 

In the Maple language (see sec. 1.4) that program-sketch becomes 

 

> H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
H := proc (K) options operator, arrow; if `and`(`and`(`or`(procname::(Not(indexed)), `<`(nops(procname), 4)), K::Ket), (op(1, K))::'identical(A, B)') then if procname::'indexed' then if nops(procname)...
 

 

Let's see it in action. Start by erasing the Physics performance remember tables, which remember results like computed before the definition of H 

 

> Library:-Forget()
 

> Typesetting:-delayDotProduct(H, Ket(A, k))
 

H[k] (95)
 

Recalling that H is Hermitian, 

> Typesetting:-delayDotProduct(Bra(B, j), H)
 

H[j] (96)
 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(B, j), H), Ket(A, k))
 

H[j, k] (97)
 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(B, j), H), Ket(A, k)), Ket(B, l))
 

H[j, k, l] (98)
 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(A, i), Bra(B, j)), H), Ket(A, k)), Ket(B, l))
 

`ℍ`[i, j, k, l] (99)
 

Note that the definition of H as a procedure does not interfere with the setting of an bracket rule for it with Ket(C, a, b), that is still working 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(C, i, j), H), Ket(C, k, l))
 

`ℋ`[i, j, k, l] (100)
 

where â„‹ = ℍ when Ket(C, k, l) is a product state. The definition of H takes precedence, so if in that definition you indicate what to do with a C Ket, that will be taken into account before the bracket rule. 

 

  • You can set the projectors for all these operators / spaces. For example,
 

> `𝕀__A` := Projector(Ket(A, i)); 1; `𝕀__B` := Projector(Ket(B, i)); 1; `𝕀__C` := Projector(Ket(C, a, b))
 

 

 

Typesetting:-mprintslash([`𝕀__A` := Sum(Physics:-`*`(Physics:-Ket(A, i), Physics:-Bra(A, i)), i = 0 .. 1)], [Sum(Physics:-`*`(Physics:-Ket(A, i), Physics:-Bra(A, i)), i = 0 .. 1)])
Typesetting:-mprintslash([`𝕀__B` := Sum(Physics:-`*`(Physics:-Ket(B, i), Physics:-Bra(B, i)), i = 0 .. 1)], [Sum(Physics:-`*`(Physics:-Ket(B, i), Physics:-Bra(B, i)), i = 0 .. 1)])
Typesetting:-mprintslash([`𝕀__C` := Sum(Sum(Physics:-`*`(Physics:-Ket(C, a, b), Physics:-Bra(C, a, b)), a = 0 .. 1), b = 0 .. 1)], [Sum(Sum(Physics:-`*`(Physics:-Ket(C, a, b), Physics:-Bra(C, a, ... (101)
 

Since the algebra rules for computing with eigenkets of A, B and C were already set in Setup(Physics:-Bracket(Physics:-Bra(A, k), Physics:-Ket(C, a, b)) = Sum(`*`(M[a, k, b, p], `*`(Physics:-Ket(B, p))), p = 0 .. 1), Physics:-Bracket(Physics:-Bra(B, k), Physics:-Ket(C, a, b)) = Sum(`*`(..., from the projectors above you can construct any subspace projector, for example 

> Typesetting:-delayDotProduct(Bra(A, m), `𝕀__C`)
 

Typesetting:-mprintslash([Sum(Sum(Sum(`*`(M[a, m, b, p], `*`(Physics:-`*`(Physics:-Ket(B, p), Physics:-Bra(C, a, b)))), p = 0 .. 1), a = 0 .. 1), b = 0 .. 1)], [Sum(Sum(Sum(`*`(M[a, m, b, p], `*`(Phys... (102)
 

> Typesetting:-delayDotProduct(`𝕀__C`, Ket(A, m))
 

Typesetting:-mprintslash([Sum(Sum(Sum(`*`(conjugate(M[a, m, b, p]), `*`(Physics:-`*`(Physics:-Ket(C, a, b), Physics:-Bra(B, p)))), p = 0 .. 1), a = 0 .. 1), b = 0 .. 1)], [Sum(Sum(Sum(`*`(conjugate(M[... (103)
 

The conjugate of M[a, m, b, p] is due to the contraction or attachment from the right of Typesetting:-delayDotProduct(Bra(A, m), `𝕀__C`), that is with 

> Dagger(Physics:-Ket(C, a, b) = Sum(Sum(`*`(M[a, j, b, p], `*`(Physics:-`*`(Physics:-Ket(A, j), Physics:-Ket(B, p)))), j = 0 .. 1), p = 0 .. 1))
 

Typesetting:-mprintslash([Physics:-Bra(C, a, b) = Sum(Sum(`*`(conjugate(M[a, j, b, p]), `*`(Physics:-`*`(Physics:-Bra(A, j), Physics:-Bra(B, p)))), j = 0 .. 1), p = 0 .. 1)], [Physics:-Bra(C, a, b) = ... (104)
 

 

The coefficients M[a, m, b, p] satisfy constraints due to the normalization of Kets of A and B. One can derive these constraints by inserting the unit operator `𝕀__C` in the identity 

> Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Bra(A, m), Bra(B, n)), `𝕀__C`), Ket(A, r)), Ket(B, s)) = Typesetting:-delayDotP...
 

Sum(Sum(`*`(conjugate(M[a, r, b, s]), `*`(M[a, m, b, n])), a = 0 .. 1), b = 0 .. 1) = `*`(Physics:-KroneckerDelta[m, r], `*`(Physics:-KroneckerDelta[n, s])) (105)
 

Transform this result into a function P  to explore the identity further 

> P := unapply(subs(Sum = sum, Sum(Sum(`*`(conjugate(M[a, r, b, s]), `*`(M[a, m, b, n])), a = 0 .. 1), b = 0 .. 1) = `*`(Physics:-KroneckerDelta[m, r], `*`(Physics:-KroneckerDelta[n, s]))), m, n, r, s)
 

Typesetting:-mprintslash([P := proc (m, n, r, s) options operator, arrow; sum(sum(`*`(conjugate(M[a, r, b, s]), `*`(M[a, m, b, n])), a = 0 .. 1), b = 0 .. 1) = `*`(Physics:-KroneckerDelta[m, r], `*`(P... (106)
 

The first and third indices refer to the quantum numbers of A, the second and fourth to B, so the right-hand sides in the following are respectively 1 and 0 

> P(1, 0, 1, 0)
 

`+`(`*`(conjugate(M[0, 1, 0, 0]), `*`(M[0, 1, 0, 0])), `*`(conjugate(M[1, 1, 0, 0]), `*`(M[1, 1, 0, 0])), `*`(conjugate(M[0, 1, 1, 0]), `*`(M[0, 1, 1, 0])), `*`(conjugate(M[1, 1, 1, 0]), `*`(M[1, 1, 1... (107)
 

> P(1, 0, 0, 0)
 

`+`(`*`(conjugate(M[0, 0, 0, 0]), `*`(M[0, 1, 0, 0])), `*`(conjugate(M[1, 0, 0, 0]), `*`(M[1, 1, 0, 0])), `*`(conjugate(M[0, 0, 1, 0]), `*`(M[0, 1, 1, 0])), `*`(conjugate(M[1, 0, 1, 0]), `*`(M[1, 1, 1... (108)
 

To get the whole system of equations satisfied by the coefficients M[a, m, b, n], use P to construct an Array with four indices running from 0..1 

> Array(`$`(0 .. 1, 4), P)
 

Typesetting:-msub(Typesetting:-mi( (109)
 

Convert the whole Array into a set of equations 

>
 

{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
{`+`(`*`(`^`(abs(M[0, 0, 0, 0]), 2)), `*`(`^`(abs(M[1, 0, 0, 0]), 2)), `*`(`^`(abs(M[0, 0, 1, 0]), 2)), `*`(`^`(abs(M[1, 0, 1, 0]), 2))) = 1, `+`(`*`(`^`(abs(M[0, 0, 0, 1]), 2)), `*`(`^`(abs(M[1, 0, 0...
(110)
 

Coherent States in Quantum Mechanics 

  • Coherent states are among the most relevant representations for the state of a quantum system. These states, which are naturally radiated from lasers, form an overcomplete basis and minimize the quantum uncertainty between position x and momentum p, so that they satisfy the Heisenberg uncertainty principle with equality and their expectation values satisfy the classical equations of motion. In this sense, coherent states represent the best possible compromise between quantities that are incompatible with regards to Heisenberg's principle. Coherent states are widely used in quantum optics, quantum information and play a key role in the description of quantum systems that can be likened to a harmonic oscillator, as for instance electromagnetic radiation, superfluids and super-conductors.
 

 

References 

[1] Cohen-Tannoudji, C.; Diu, B.; and Laloe, F. Quantum Mechanics. Paris, France: Hermann, 1977. 

[2] Massachusetts Institute of Technology OpenCourseWare, Quantum Physics II, Quantum Dynamics. 

Definition and the basics 

> restart; 1; with(Physics); -1
 

 

Set a quantum operator A and corresponding annihilation / creation operators 

> Setup(quantumoperators = A)
 

[quantumoperators = {A}] (111)
 

> am := Annihilation(A)
 

Typesetting:-mprintslash([am := `#msup(mi( (112)
 

> ap := Creation(A)
 

Typesetting:-mprintslash([ap := `#msup(mi( (113)
 

In what follows, on the left-hand sides the product operator used is `*`, which properly represents, but does not perform the attachment of Bras, Kets, and operators. On the right-hand sides the product operator is `.`, that performs the attachments. Since the introduction of Physics in the Maple system, we have that 

> `*`(am, `*`(Ket(A, n))) = Typesetting:-delayDotProduct(am, Ket(A, n))
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (114)
 

> (%Bracket = Bracket)(Bra(A, n), Ket(A, n))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (115)
 

> (%Bracket = Bracket)(Bra(A, n), Ket(A, m))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (116)
 

New in Maple 2019: coherent states, the eigenstates of the annihilation operator `#msup(mi(, with all of their properties, are now understood as such by the system 

> `*`(am, `*`(Ket(am, alpha))) = Typesetting:-delayDotProduct(am, Ket(am, alpha))
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (117)
 

Ket(`#msup(mi( is an eigenket of `#msup(mi( but not of  `#msup(mi( 

> Typesetting:-delayDotProduct(ap, Ket(am, alpha))
 

Typesetting:-mprintslash([Physics:-`.`(`#msup(mi( (118)
 

The norm of these states is equal to 1 

> (%Bracket = Bracket)(Bra(am, alpha), Ket(am, alpha))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (119)
 

These states, however, are not orthonormal as the occupation number states Ket(A, n) are, and since `#msup(mi( is not Hermitian, its eigenvalues are not real but complex numbers. Instead of (%Bracket = Bracket)(Bra(A, n), Ket(A, m)), in Maple 2019 we have 

> (%Bracket = Bracket)(Bra(am, alpha), Ket(am, beta))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (120)
 

At alpha = beta,  

> simplify(eval(%Bracket(Physics:-Bra(`#msup(mi(
 

1 = 1 (121)
 

Their scalar product with the occupation number states Ket(A, m), using the inert %Bracket on the left-hand side and the active Bracket on the other side: 

> (%Bracket = Bracket)(Bra(A, n), Ket(am, alpha))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (122)
 

The expansion of coherent states into occupation number states, first representing the product operation using `*`, then performing the attachments replacing `*` by `.` 

> Projector(Ket(A, n), dimension = infinity)
 

Typesetting:-mprintslash([Sum(Physics:-`*`(Physics:-Ket(A, n), Physics:-Bra(A, n)), n = 0 .. infinity)], [Sum(Physics:-`*`(Physics:-Ket(A, n), Physics:-Bra(A, n)), n = 0 .. infinity)]) (123)
 

> Ket(am, alpha) = `*`(Sum(Physics:-`*`(Physics:-Ket(A, n), Physics:-Bra(A, n)), n = 0 .. infinity), `*`(Ket(am, alpha)))
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (124)
 

> eval(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (125)
 

Taking all into account, 

> Typesetting:-delayDotProduct(Bra(A, n), Physics:-Ket(`#msup(mi(
 

`/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(`^`(alpha, n))), `*`(`^`(factorial(n), `/`(1, 2)))) = `/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(`^`(alpha,... (126)
 

Hide now the ket label. When in doubt, input show to see the Kets with their labels explicitly shown 

> Setup(hide = true)
 

 

 

`*`(`* Partial match of  '`, `*`(hide, `*`(`' against keyword '`, `*`(hideketlabel, `*`(`' `)))))
_______________________________________________________
[hideketlabel = true] (127)
 

Define eigenkets of the annihilation operator, with two different eigenvalues for experimentation 

> `K__α` := Ket(am, alpha)
 

Typesetting:-mprintslash([`K__α` := Physics:-Ket(`#msup(mi( (128)
 

> `K__β` := Ket(am, beta)
 

Typesetting:-mprintslash([`K__β` := Physics:-Ket(`#msup(mi( (129)
 

Because the properties of coherent states are now known to the system, the following computations proceed automatically in Maple 2019. The left-hand sides use the `*`, while the right-hand sides use the `.`  

> (`*` = `.`)(Dagger(`K__α`), ap, am, `K__α`)
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (130)
 

> (`*` = `.`)(Dagger(`K__α`), `+`(ap, am), `K__α`)
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (131)
 

> (`*` = `.`)(Dagger(`K__α`), `+`(ap, `-`(am)), `K__α`)
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (132)
 

> (`*` = `.`)(Dagger(`K__α`), `*`(`^`(`+`(ap, am), 2)), `K__α`)
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (133)
 

Properties of Coherent states 

The mean value of the occupation number N 

 

The occupation number operator N is given by 

> N := Typesetting:-delayDotProduct(ap, am)
 

Typesetting:-mprintslash([N := Physics:-`*`(`#msup(mi( (134)
 

`*`(N, `*`(` is Hermitian`)) 

> %Dagger(N) = N
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (135)
 

> value(%Dagger(Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (136)
 

N is diagonal in the Ket(A, n) basis of the Fock (occupation number) space 

> (`*` = `.`)(Bra(A, n), N, Ket(A, p))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(A, n), `#msup(mi( (137)
 

  • The mean value of N in a coherent state `≡`(Ket(`#msup(mi(
 

> Bracket(%N)[alpha] = %Bracket(Bra(am, alpha), N, Ket(am, alpha))
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mrow(Typesetting:-mi( (138)
 

> value(Physics:-Bracket(%N)[alpha] = %Bracket(Physics:-Bra(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Bracket(%N)[alpha] = `*`(`^`(abs(alpha), 2))], [Physics:-Bracket(%N)[alpha] = `*`(`^`(abs(alpha), 2))]) (139)
 

The mean value of `*`(`^`(N, 2)) 

> Bracket(`*`(`^`(%N, 2)))[alpha] = %Bracket(Bra(am, alpha), `*`(`^`(N, 2)), Ket(am, alpha))
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mrow(Typesetting:-mi( (140)
 

> value(Physics:-Bracket(`*`(`^`(%N, 2)))[alpha] = %Bracket(Physics:-Bra(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Bracket(`*`(`^`(%N, 2)))[alpha] = `+`(`*`(`^`(abs(alpha), 4)), `*`(`^`(abs(alpha), 2)))], [Physics:-Bracket(`*`(`^`(%N, 2)))[alpha] = `+`(`*`(`^`(abs(alpha), 4)), `*... (141)
 

The standard deviation `ΔN` = sqrt(`+`(`-`(`*`(`^`(Bracket(%N)[alpha], 2))), Bracket(`^`(%N, 2))[alpha])) for a state Ket(`#msup(mi( 

> `*`(`^`(`+`(Physics:-Bracket(`*`(`^`(%N, 2)))[alpha] = `+`(`*`(`^`(abs(alpha), 4)), `*`(`^`(abs(alpha), 2))), `-`(`^`(Physics:-Bracket(%N)[alpha] = `*`(`^`(abs(alpha), 2)), 2))), `/`(1, 2)))
 

`*`(`^`(`+`(`-`(`*`(`^`(Physics:-Bracket(%N)[alpha], 2))), Physics:-Bracket(`*`(`^`(%N, 2)))[alpha]), `/`(1, 2))) = abs(alpha) (142)
 

In conclusion, a coherent state has a finite spreading `ΔN` = abs(alpha).  Coherent states are good approximations for the states of a laser, where the laser intensity I  is proportional to the mean value of the photon number, I ∝ Bracket(%N)[alpha] = `*`(`^`(abs(alpha), 2)), and so the intensity fluctuation, `∝`(sqrt(I), abs(alpha)). 

  • The mean value of the occupation number N in an occupation number state `≡`(Ket(A, n), Ket(n))
 

> Bracket(%N)[n] = %Bracket(Bra(A, n), N, Ket(A, n))
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mrow(Typesetting:-mi( (143)
 

> value(Physics:-Bracket(%N)[n] = %Bracket(Physics:-Bra(A, n), Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Bracket(%N)[n] = n], [Physics:-Bracket(%N)[n] = n]) (144)
 

The mean value of the occupation number N in a state Ket(A, n) is thus n itself, as expected since Ket(A, n)represents a (Fock space) state of n (quasi-) particles. Accordingly, 

> Bracket(`*`(`^`(%N, 2)))[n] = %Bracket(Bra(A, n), `*`(`^`(N, 2)), Ket(A, n))
 

(145)
 

> value(Physics:-Bracket(`*`(`^`(%N, 2)))[n] = %Bracket(Physics:-Bra(A, n), Physics:-`^`(Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Bracket(`*`(`^`(%N, 2)))[n] = `*`(`^`(n, 2))], [Physics:-Bracket(`*`(`^`(%N, 2)))[n] = `*`(`^`(n, 2))]) (146)
 

The standard deviation `ΔN` = sqrt(`+`(`-`(`*`(`^`(Bracket(%N)[n], 2))), Bracket(`^`(%N, 2))[n])) for a state Ket(`#msup(mi(, is thus 

> `*`(`^`(`+`(Physics:-Bracket(`*`(`^`(%N, 2)))[n] = `*`(`^`(n, 2)), `-`(`^`(Physics:-Bracket(%N)[n] = n, 2))), `/`(1, 2)))
 

`*`(`^`(`+`(`-`(`*`(`^`(Physics:-Bracket(%N)[n], 2))), Physics:-Bracket(`*`(`^`(%N, 2)))[n]), `/`(1, 2))) = 0 (147)
 

That is, in a Fock state, `ΔN` = 0, there is no intensity fluctuation. 

 

 

The specific properties of coherent states implemented in Maple 2019 can be derived explicitly departing from the projection of Ket(`#msup(mi( into the Ket(A, m)basis of occupation number states and the definition of `#msup(mi( as the operator that annihilates the vacuum `#msup(mi(Ket(A, n) = 0 

> Ket(am, alpha) = `*`(Sum(Physics:-`*`(Physics:-Ket(A, n), Physics:-Bra(A, n)), n = 0 .. infinity), `*`(Ket(am, alpha)))
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (148)
 

> eval(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (149)
 

To derive `*`(`#msup(mi( from the formula above, start multiplying by `#msup(mi( 

> `*`(am, `*`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (150)
 

In view of `*`(`#msup(mi(, discard the first term of the sum 

> subs(0 = 1, Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (151)
 

Change variables n = `+`(k, 1); in the result rename proc (k) options operator, arrow; n end proc 

> subs(k = n, PDEtools:-dchange(n = `+`(k, 1), Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (152)
 

Activate the product `*`(`#msup(mi( by replacing, in the right-hand side, the product operator `*` by `.` 

> lhs(Physics:-`*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`#msup(mi( (153)
 

By inspection the right-hand side of lhs(Physics:-`*`(`#msup(mi( is equal to alpha times the right-hand side of eval(Physics:-Ket(`#msup(mi( 

> `+`(`*`(alpha, `*`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([`+`(`*`(alpha, `*`(Physics:-Ket(`#msup(mi( (154)
 

> combine(`+`(`*`(alpha, `*`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([`+`(`*`(alpha, `*`(Physics:-Ket(`#msup(mi( (155)
 

  • Overview of the coherent states distribution
 

 

Consider the projection of Ket(`#msup(mi( over an occupation number state Ket(A, n) 

> %Bracket(Bra(A, n), lhs(Physics:-Ket(`#msup(mi(
 

(156)
 

An overview of the distribution of coherent states Ket(`#msup(mi( for a sample of values of n and alpha is thus as follows 

> plot3d(rhs(%Bracket(Physics:-Bra(A, n), Physics:-Ket(`#msup(mi(
 

Plot_2d
 

The distribution can be explored for ranges of values of n and alpha using Explore 

> NA := Typesetting:-Typeset(Bracket(Bra(A, n), Ket(am, alpha))); -1
 

> Explore(plot(rhs(%Bracket(Physics:-Bra(A, n), Physics:-Ket(`#msup(mi(
 

Embedded component

alpha Embedded component

5

 

 

To verify this identity, construct each of the three terms, then simplify the result. Recalling the projection of a coherent stateKet(`#msup(mi(into the Ket(A, m)basis of occupation number states, 

> Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (157)
 

The third term of this identity is thus 

> T__3 := `+`(`*`(`/`(1, 2), `*`(conjugate(alpha), `*`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([T__3 := `+`(`*`(`/`(1, 2), `*`(conjugate(alpha), `*`(Physics:-Ket(`#msup(mi( (158)
 

The first term, on the left-hand side, is obtained multiplying by `#msup(mi( 

> Typesetting:-delayDotProduct(ap, Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`#msup(mi( (159)
 

To have the three terms with Ket(A, n)in the summand, change variables n = `+`(k, `-`(1)); in the result rename proc (k) options operator, arrow; n end proc 

> subs(k = n, convert(PDEtools:-dchange(n = `+`(k, `-`(1)), Physics:-`.`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`#msup(mi( (160)
 

The radical in the summand can be rewritten taking into account that, when n is a positive integer, 

> `/`(1, `*`(sqrt(factorial(`+`(n, `-`(1)))))) = `/`(`*`(sqrt(n)), `*`(sqrt(factorial(n))))
 

`/`(1, `*`(`^`(factorial(`+`(n, `-`(1))), `/`(1, 2)))) = `/`(`*`(`^`(n, `/`(1, 2))), `*`(`^`(factorial(n), `/`(1, 2)))) (161)
 

This identity can be verified as follows 

> `assuming`([simplify((`+`(lhs, `-`(rhs)))(`/`(1, `*`(`^`(factorial(`+`(n, `-`(1))), `/`(1, 2)))) = `/`(`*`(`^`(n, `/`(1, 2))), `*`(`^`(factorial(n), `/`(1, 2))))))], [n::posint])
 

0 (162)
 

> subs(`/`(1, `*`(`^`(factorial(`+`(n, `-`(1))), `/`(1, 2)))) = `/`(`*`(`^`(n, `/`(1, 2))), `*`(`^`(factorial(n), `/`(1, 2)))), Physics:-`.`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`#msup(mi( (163)
 

The summand, at n = 0, is equal to 0 

> summand := op(1, rhs(Physics:-`.`(`#msup(mi(
 

Typesetting:-mprintslash([summand := `/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(`^`(alpha, `+`(n, `-`(1))), `*`(n, `*`(Physics:-Ket(A, n))))), `*`(`^`(factorial(n), `/`(1, 2))... (164)
 

> eval(summand, n = 0)
 

0 (165)
 

Rewriting then the sum to start from 0 

> T__1 := convert(subs(1 .. infinity = 0 .. infinity, Physics:-`.`(`#msup(mi(
 

Typesetting:-mprintslash([T__1 := Physics:-`.`(`#msup(mi( (166)
 

 

The second term of this identity is obtained by differentiation of Physics:-Ket(`#msup(mi( 

> T__2 := diff(Physics:-Ket(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (167)
 

Putting the three terms together, 

> `+`(T__1, `-`(T__2), `-`(T__3))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(168)
 

Combining the sums the identity is verified 

> combine(`+`(Physics:-`.`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (169)
 

`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi( =  

The coherent state can be constructed from the vacuum state using the operator `𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi(. 

> Setup(quantumoperators = {`𝒟`})
 

[quantumoperators = {`𝒟`, A}] (170)
 

> `𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(ap)))))
 

`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi( (171)
 

New in Maple 2019, the conversion network for mathematical functions can be used with not-commutative variables;  develop the exponential function of `#msup(mi( in power series: 

> exp(`*`(alpha, `*`(ap))) = convert(exp(`*`(alpha, `*`(ap))), Sum)
 

exp(`*`(alpha, `*`(`#msup(mi( (172)
 

So `𝒟`[1]becomes 

> subs(exp(`*`(alpha, `*`(`#msup(mi(
 

`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(Sum(`/`(`*`(Physics:-`^`(`*`(alpha, `*`(`#msup(mi( (173)
 

Therefore, for `𝒟`[1]·, we have 

> Typesetting:-delayDotProduct(`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(Sum(`/`(`*`(Physics:-`^`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Sum(`/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(Physics:-`.`(Physics:-`^`(`*`(alpha, `*`(`#msup(mi( (174)
 

> combine(Expand(Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Sum(`/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(Physics:-`.`(Physics:-`^`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Sum(`/`(`*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(`^`(alpha, n), `*`(Physics:-Ket(A, n)))), `*`(`^`(fact... (175)
 

By inspection the right-hand side is already the projection of Ket(`#msup(mi(into the basis or occupation number states Ket(A, n)computed previously in eval(Physics:-Ket(`#msup(mi( 

> Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (176)
 

> `+`(`*`(`𝒟`[1], `*`(Physics:-Ket(A, 0))), `-`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([`+`(Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)), `-`(Physics:-Ket(`#msup(mi( (177)
 

> isolate(`+`(Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)), `-`(Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Physics:-Ket(`#msup(mi( (178)
 

Remark: `𝒟`[1] is not unitary, `<>`(`*`(`𝒟`[1], `*`(`#msup(mi( 

> `*`(Dagger(`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (179)
 

> simplify(Physics:-`*`(Physics:-Dagger(`𝒟`[1]), `𝒟`[1]) = `*`(`^`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), 2), `*`(Physics:-`*`(exp(`*`(conjugate(alpha), `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (180)
 

exp(`+`(`*`(alpha, `*`(`#msup(mi( =  

 

Here, we use another operator, `𝒟`[2]= exp(`+`(`*`(alpha, `*`(`#msup(mi( to construct from the vacuum. `𝒟`[2] is sometime called the "displacement" operator. It has the advantage over `𝒟`[1] that `𝒟`[2]is unitary,`and`(`*`(`𝒟`[2], `*`(`#msup(mi(As a consequence:  . 

> `𝒟`[2] = exp(`+`(`*`(alpha, `*`(ap)), `-`(`*`(conjugate(alpha), `*`(am)))))
 

`𝒟`[2] = exp(`+`(`*`(alpha, `*`(`#msup(mi( (181)
 

This operator is unitary 

> `*`(Dagger(`𝒟`[2] = exp(`+`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (182)
 

> Simplify(Physics:-`*`(`𝒟`[2], Physics:-Dagger(`𝒟`[2])) = Physics:-`*`(exp(`+`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (183)
 

To verify that `*`(`𝒟`[2], `*`(`#msup(mi(one can proceed as in the above, or directly compute their commutator, expecting Physics:-Commutator(`𝒟`[2], `#msup(mi(. 

> %Commutator(`𝒟`[2] = exp(`+`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mi( (184)
 

> value(%)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (185)
 

For `𝒟`[2]·, start multiplying by `𝒟`[2] 

> `*`(Ket(A, 0), `*`(`𝒟`[2])) = `*`(Ket(A, 0), `*`(exp(`+`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`*`(exp(`+`(`*`(alpha, `*`(`#msup(mi( (186)
 

> simplify(expand(Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`*`(exp(`+`(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = `*`(Physics:-`*`(exp(`*`(alpha, `*`(`#msup(mi( (187)
 

Recalling the definition of `𝒟`[1]in the previous section 

> `𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi(
 

`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi( (188)
 

The expression above simplify(expand(Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`*`(exp(`+`(`*`(alpha, `*`(`#msup(mi( can be simplified using this definition 

> isolate(`𝒟`[1] = `*`(exp(`+`(`-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))), `*`(exp(`*`(alpha, `*`(`#msup(mi(
 

exp(`*`(alpha, `*`(`#msup(mi( (189)
 

> simplify(subs(exp(`*`(alpha, `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`*`(`𝒟`[1], exp(`+`(`-`(`*`(conjugate(alpha), `*`(`#msup(mi( (190)
 

In turn `*`(exp(`+`(`-`(`*`(conjugate(alpha), `*`(`#msup(mi( can be computed by replacing the `*` product by a dot product `.` 

> eval(Physics:-`*`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`*`(`𝒟`[1], exp(`+`(`-`(`*`(conjugate(alpha), `*`(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0))], [Physics:-`.`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-`.`(`𝒟`[1], Phy... (191)
 

Finally, we arrive at the desired result recalling the result of the previous section, 

> Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Physics:-Ket(`#msup(mi( (192)
 

> subs(Physics:-`.`(`𝒟`[1], Physics:-Ket(A, 0)) = Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`.`(`𝒟`[2], Physics:-Ket(A, 0)) = Physics:-Ket(`#msup(mi( (193)
 

`<|>`(beta, alpha) = exp(`+`(`*`(conjugate(beta), `*`(alpha)), `-`(`*`(`/`(1, 2), `*`(`^`(abs(beta), 2)))), `-`(`*`(`/`(1, 2), `*`(`^`(abs(alpha), 2)))))) 

 

The identity in the title can be derived departing again from the projection of a coherent stateKet(`#msup(mi(into the Ket(A, m)basis of occupation number states 

> Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Ket(`#msup(mi( (194)
 

> Dagger(subs({alpha = beta, n = k}, Physics:-Ket(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-Bra(`#msup(mi( (195)
 

Taking the `*` product of these two expressions 

> `*`(Physics:-Bra(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (196)
 

Perform the attachment of Bras and Kets on the right-hand side by replacing `*` by `.`, evaluating the sum and simplifying the result 

> lhs(Physics:-`*`(Physics:-Bra(`#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Bra(`#msup(mi( (197)
 

  • Overview of the real and imaginary part of `<|>`(beta, alpha)
 

 

In most cases, alpha and beta are complex valued numbers. Below, the plots assume that alpha and beta are both real. To take into account the general case, the possibility to tune a phase difference theta between alpha and beta is explicitly added, so that lhs(Physics:-`*`(Physics:-Bra(`#msup(mi( becomes 

 

> %Bracket(Bra(am, beta), Ket(am, alpha)) = subs(conjugate(beta) = `*`(conjugate(beta), `*`(exp(`*`(I, `*`(theta))))), rhs(Physics:-`*`(Physics:-Bra(`#msup(mi(
 

(198)
 

> Explore(plot3d(Re(rhs(%Bracket(Physics:-Bra(`#msup(mi(
Explore(plot3d(Re(rhs(%Bracket(Physics:-Bra(`#msup(mi(
 

Embedded component

theta Embedded component

0.314

 

The Zassenhaus formula and the algebra of the Pauli matrices 

 

The implementation of the Pauli matrices and their algebra were reviewed for Maple 2019, including the algebraic manipulation of nested commutators, resulting in faster computations using simpler and more flexible input. As it frequently happens, improvements of this type suddenly transform research problems presented in the literature as untractable in practice, into tractable. 

As an illustration, we tackle below the derivation of the coefficients entering the Zassenhaus formula shown in section 4 of [1] for the Pauli matrices up to order 10 (results in the literature go up to order 5). The computation presented can be reused to compute these coefficients up to any desired higher-order (hardware limitations may apply). A number of examples which exploit this formula and its dual, the Baker-Campbell-Hausdorff formula, occur in connection with the Weyl prescription for converting a classical function to a quantum operator (see sec. 5 of [1]), as well as when solving the eigenvalue problem for classes of mathematical-physics partial differential equations [2].  

References 

  • [1] R.M. Wilcox. "Exponential Operators and Parameter Differentiation in Quantum Physics", Journal of Mathematical Physics, V.8, 4, (1967.
 

  • [2] S. Steinberg. "Applications of the lie algebraic formulas of Baker, Campbell, Hausdorff, and Zassenhaus to the calculation of explicit solutions of partial differential equations". Journal of Differential Equations, V.26, 3, 1977.
 

  • [3] K. Huang. "Statistical Mechanics". John Wiley & Sons, Inc. 1963, p217, Eq.(10.60).
 

 

Formulation of the problem 

The Zassenhaus formula expresses exp(`*`(lambda, `*`(`+`(A, B)))) as an infinite product of exponential operators involving nested commutators of increasing complexity 


 

Given A, B and their commutator E = %Commutator(A, B), if A and B commute with E, C[n] = 0 for `>=`(n, 3) and the Zassenhaus formula reduces to the product of the first three exponentials above. The interest here is in the general case, when `<>`(%Commutator(A, E), 0) and `<>`(%Commutator(B, E), 0), and the goal is to compute the Zassenhaus coefficients C[n]in terms of A, B for arbitrary n. Following [1], in that general case, differentiating the Zassenhaus formula with respect to lambda and multiplying from the right by exp(`+`(`-`(`*`(lambda, `*`(`+`(A, B)))))) one obtains 

 

This is an intricate formula, which however (see eq.(4.20) of [1]) can be represented in abstract form as 

 

 

from where an equation to be solved for each C[n] is obtained by equating to 0 the coefficient of `^`(lambda, `+`(n, `-`(1))). In this formula, the repeated commutator bracket is defined inductively in terms of the standard commutator %Commutator(A, B)by 

{B, `^`(A, 0)} = B, {B, `^`(A, `+`(n, 1))} = %Commutator(A, {`^`(A, n), `^`(B, n)}) 

{C[j], `^`(B, n), `^`(A, 0)} = {C[j], `^`(B, n)}, {C[j], `^`(A, m), `^`(B, n)} = %Commutator(A, {`^`(A, `-`(m, 1)), `^`(B, n), `^`(C[j], k)}) 

and higher-order repeated-commutator brackets are similarly defined. For example, taking the coefficient of lambda and `*`(`^`(lambda, 2)) and respectively solving each of them for C[2] and C[3] one obtains  

C[2] = `+`(`-`(`*`(`/`(1, 2), `*`(%Commutator(A, B))))) 

C[3] = `+`(`*`(`/`(1, 2), `*`(%Commutator(B, %Commutator(A, B))))) 

This method is used in [3] to treat quantum deviations from the classical limit of the partition function for both a Bose-Einstein and Fermi-Dirac gas. The complexity of the computation of C[n] grows rapidly and in the literature only the coefficients up to C[5] have been published. Taking advantage of developments in the Physics package for Maple 2019, below we show the computation up to C[10] and provide a compact approach to compute them up to arbitrary order. 

 

Computing up to C[10] 

Set the signature of spacetime such that its space part is equal to +++ and use lowercaselatin letters to represent space indices. Set also A, B and C[n] to represent quantum operators 

> restart; 1; with(Physics); -1
 

> Setup(op = {A, B, C}, signature = `+++-`, spaceindices = lowercaselatin)
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[quantumoperators = {A, B, C}, signature = `+ + + -`, spaceindices = lowercaselatin] (199)
 

To illustrate the computation up to C[10], a convenient example, where the commutator algebra is closed, consists of taking A and B as Pauli Matrices which, multiplied by the imaginary unit, form a basis for the `𝔰𝔲`(2)group, which in turn exponentiate to the relevant Special Unitary Group SU(2). The algebra for the Pauli matrices involves a commutator and an anticommutator 

> Library:-DefaultAlgebraRules(Psigma)
 

Typesetting:-mprintslash([Typesetting:-mcomplete(Typesetting:-msub(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (200)
 

Assign now A and B to two Pauli matrices, for instance 

> A := Psigma[1]
 

Typesetting:-mprintslash([A := Physics:-Psigma[1]], [Physics:-Psigma[1]]) (201)
 

> B := Psigma[3]
 

Typesetting:-mprintslash([B := Physics:-Psigma[3]], [Physics:-Psigma[3]]) (202)
 

Next, to extract the coefficient of `^`(lambda, n) from 

 

to solve it for C[`+`(n, 1)] we note that each term has a factor `^`(lambda, m) multiplying a sum, so we only need to take into account the first `+`(n, 1) terms (sums) and in each sum replace infinity by the corresponding `+`(n, `-`(m)). For example, given to compute C[3] we only need to compute these first three terms: 

0 = `+`(Sum(`/`(`*`(`^`(lambda, n), `*`({B, `^`(A, n)})), `*`(factorial(n))), n = 1 .. 2), `*`(2, `*`(lambda, `*`(Sum(Sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`({C[2], `^`(A, m), `^`(B, n)})), `*`(factor... 

then solving for C[3] one gets C[3] = `+`(`*`(`/`(1, 3), `*`(%Commutator(B, %Commutator(A, B)))), `*`(`/`(1, 6), `*`(%Commutator(A, %Commutator(A, B))))). 

Also, since to compute C[n] we only need the coefficient of `^`(lambda, `+`(n, `-`(1))), it is not necessary to compute all the terms of each multiple-sum. One way of restricting the multiple-sums to only one power of lambda consists of using multi-index summation, available in the Physics package. For that purpose, redefine sum to extend its functionality with multi-index summation 

> Setup(redefinesum = true)
 

[redefinesum = true] (203)
 

Now we can represent the same computation of C[3] without multiple sums and without computing unnecessary terms as  

0 = `+`(Sum(`/`(`*`(`^`(lambda, n), `*`({B, `^`(A, n)})), `*`(factorial(n))), n = 1), `*`(2, `*`(lambda, `*`(Sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`({C[2], `^`(A, m), `^`(B, n)})), `*`(factorial(n), `... 

Finally, we need a computational representation for the repeated commutator bracket  

{B, `^`(A, 0)} = B, {B, `^`(A, `+`(n, 1))} = %Commutator(A, {`^`(A, n), `^`(B, n)}) 

One way of representing this commutator bracket operation is defining a procedure, say F, with a cache to avoid recomputing lower order nested commutators, as follows 

> F := proc (A, B, n) option cache; if n::negint then 0 elif n = 0 then B elif n::posint then %Commutator(A, F(A, B, `+`(n, `-`(1)))) else 'F(A, B, n)' end if end proc; 1
F := proc (A, B, n) option cache; if n::negint then 0 elif n = 0 then B elif n::posint then %Commutator(A, F(A, B, `+`(n, `-`(1)))) else 'F(A, B, n)' end if end proc; 1
 

Typesetting:-mprintslash([F := proc (A, B, n) option cache; if n::negint then 0 elif n = 0 then B elif n::posint then %Commutator(A, F(A, B, `+`(n, `-`(1)))) else 'F(A, B, n)' end if end proc], [proc ...
Typesetting:-mprintslash([F := proc (A, B, n) option cache; if n::negint then 0 elif n = 0 then B elif n::posint then %Commutator(A, F(A, B, `+`(n, `-`(1)))) else 'F(A, B, n)' end if end proc], [proc ...
Typesetting:-mprintslash([F := proc (A, B, n) option cache; if n::negint then 0 elif n = 0 then B elif n::posint then %Commutator(A, F(A, B, `+`(n, `-`(1)))) else 'F(A, B, n)' end if end proc], [proc ...
(204)
 

For example, 

> F(A, B, 1)
 

Typesetting:-mrow(Typesetting:-mi( (205)
 

> F(A, B, 2)
 

Typesetting:-mrow(Typesetting:-mi( (206)
 

> F(A, B, 3)
 

Typesetting:-mrow(Typesetting:-mi( (207)
 

We can set now the value of C[2] 

> C[2] := `+`(`-`(`*`(`/`(1, 2), `*`(Commutator(A, B)))))
 

Typesetting:-mprintslash([C[2] := `*`(I, `*`(Physics:-Psigma[2]))], [`*`(I, `*`(Physics:-Psigma[2]))]) (208)
 

and enter the formula that involves only multi-index summation 

> H := `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = 2), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*`(factor...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn( (209)
 

from where we compute C[3] by solving for it the coefficient of `*`(`^`(lambda, 2)), and since due to the multi-index summation this expression already contains `*`(`^`(lambda, 2)) as a factor,  

> C[3] = solve(value(H), C[3])
 

C[3] = `+`(`-`(`*`(`/`(4, 3), `*`(Physics:-Psigma[1]))), `*`(`/`(2, 3), `*`(Physics:-Psigma[3]))) (210)
 

In order to generalize the formula for H for higher powers of lambda, the right-hand side of the multi-index summation limit can be expressed in terms of an abstract N, and H transformed into a mapping: 

 

> H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
 

Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(`/`(`*`(`^`(lambda, N), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], N))), `*`(factorial(N))), `*`(2, `*`(lambda, `*`(sum(Physics:...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(`/`(`*`(`^`(lambda, N), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], N))), `*`(factorial(N))), `*`(2, `*`(lambda, `*`(sum(Physics:...
(211)
 

Now we have 

> H(0)
 

Physics:-Psigma[3] (212)
 

> H(1)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (213)
 

The following is already equal to H := `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = 2), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*`(factor... 

> H(2)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn( (214)
 

In this way, we can reproduce the results published in the literature for the coefficients of Zassenhaus formula up to C[5] by adding two more multi-index sums to H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*.... Unassign C first
 

> unassign(C)
 

> H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
H := unapply(`+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(A, B, n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(sum(`/`(`*`(`^`(lambda, `+`(n, m)), `*`(F(A, F(B, C[2], n), m))), `*`(factorial(n), `*...
 

We compute now up to C[5] in one go 

> for j to 4 do C[`+`(j, 1)] := solve(value(H(j)), C[`+`(j, 1)]) end do; 1
 

 

 

 

Typesetting:-mprintslash([C[2] := `*`(I, `*`(Physics:-Psigma[2]))], [`*`(I, `*`(Physics:-Psigma[2]))])
Typesetting:-mprintslash([C[3] := `+`(`-`(`*`(`/`(4, 3), `*`(Physics:-Psigma[1]))), `*`(`/`(2, 3), `*`(Physics:-Psigma[3])))], [`+`(`-`(`*`(`/`(4, 3), `*`(Physics:-Psigma[1]))), `*`(`/`(2, 3), `*`(Phy...
Typesetting:-mprintslash([C[4] := `+`(`*`(`*`(`/`(4, 3), `*`(I)), `*`(Physics:-Psigma[2])), Physics:-Psigma[1], `*`(2, `*`(Physics:-Psigma[3])))], [`+`(`*`(`*`(`/`(4, 3), `*`(I)), `*`(Physics:-Psigma[...
Typesetting:-mprintslash([C[5] := `+`(`-`(`*`(`+`(`*`(`/`(16, 3), `*`(I))), `*`(Physics:-Psigma[2]))), `-`(`*`(`/`(8, 9), `*`(Physics:-Psigma[1]))), `-`(`*`(`/`(158, 45), `*`(Physics:-Psigma[3]))))], ... (215)
 

The nested-commutator expression solved in the last step for C[5] is 

> H(4)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfrac(Typesetting:-mn(
(216)
 

With everything understood, we want now to extend these results generalizing them into an approach to compute an arbitrarily large coefficient C[n], then use that generalization to compute all the Zassenhaus coefficients up to C[10]. To type the formula for H for higher powers of lambda is however prone to typographical mistakes. The following is a program, using the Maple programming language, that produces these formulas for an arbitrary integer power of lambda: 

Embedded component

 

This Formula program uses a sequence of summation indices with as much indices as the order of the coefficient C[n] we want to compute, in this case we need 10 of them 

> summation_indices := n, m, k, l, p, q, r, s, t, u
 

Typesetting:-mprintslash([summation_indices := n, m, k, l, p, q, r, s, t, u], [n, m, k, l, p, q, r, s, t, u]) (217)
 

To avoid interference of the results computed in the loop for j to 4 do C[`+`(j, 1)] := solve(value(H(j)), C[`+`(j, 1)]) end do; 1, unassign C again 

> unassign(C)
 

 

Now the formulas typed by hand, used lines above to compute each of C[2], C[3] and C[5], are respectively constructed by the computer 

> Formula(A, B, C, 2)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (218)
 

> Formula(A, B, C, 3)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (219)
 

> Formula(A, B, C, 5)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi(
(220)
 

 

Construct then the formula for C[10] and make it be a mapping with respect to N, as done for C[5] after H(2) 

> H := unapply(Formula(A, B, C, 10), N)
 

Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
Typesetting:-mprintslash([H := proc (N) options operator, arrow; `+`(sum(`/`(`*`(`^`(lambda, n), `*`(F(Physics:-Psigma[1], Physics:-Psigma[3], n))), `*`(factorial(n))), n = N), `*`(2, `*`(lambda, `*`(...
(221)
 

Compute now the coefficients of the Zassenhaus formula up to C[10] all in one go 

> for j to 9 do C[`+`(j, 1)] := solve(value(H(j)), C[`+`(j, 1)]) end do; 1
 

 

 

 

 

 

 

 

 

Typesetting:-mprintslash([C[2] := `*`(I, `*`(Physics:-Psigma[2]))], [`*`(I, `*`(Physics:-Psigma[2]))])
Typesetting:-mprintslash([C[3] := `+`(`-`(`*`(`/`(4, 3), `*`(Physics:-Psigma[1]))), `*`(`/`(2, 3), `*`(Physics:-Psigma[3])))], [`+`(`-`(`*`(`/`(4, 3), `*`(Physics:-Psigma[1]))), `*`(`/`(2, 3), `*`(Phy...
Typesetting:-mprintslash([C[4] := `+`(`*`(`*`(`/`(4, 3), `*`(I)), `*`(Physics:-Psigma[2])), Physics:-Psigma[1], `*`(2, `*`(Physics:-Psigma[3])))], [`+`(`*`(`*`(`/`(4, 3), `*`(I)), `*`(Physics:-Psigma[...
Typesetting:-mprintslash([C[5] := `+`(`-`(`*`(`+`(`*`(`/`(16, 3), `*`(I))), `*`(Physics:-Psigma[2]))), `-`(`*`(`/`(8, 9), `*`(Physics:-Psigma[1]))), `-`(`*`(`/`(158, 45), `*`(Physics:-Psigma[3]))))], ...
Typesetting:-mprintslash([C[6] := `+`(`*`(`*`(`/`(1078, 405), `*`(I)), `*`(Physics:-Psigma[2])), `*`(`/`(1030, 81), `*`(Physics:-Psigma[1])), `-`(`*`(`/`(8, 81), `*`(Physics:-Psigma[3]))))], [`+`(`*`(...
Typesetting:-mprintslash([C[7] := `+`(`*`(`*`(`/`(11792, 243), `*`(I)), `*`(Physics:-Psigma[2])), `*`(`/`(358576, 42525), `*`(Physics:-Psigma[1])), `*`(`/`(12952, 135), `*`(Physics:-Psigma[3])))], [`+...
Typesetting:-mprintslash([C[8] := `+`(`*`(`*`(`/`(35837299048, 17222625), `*`(I)), `*`(Physics:-Psigma[2])), `*`(`/`(87277417, 492075), `*`(Physics:-Psigma[1])), `*`(`/`(833718196, 820125), `*`(Physic...
Typesetting:-mprintslash([C[9] := `+`(`-`(`*`(`+`(`*`(`/`(449018539801088, 104627446875), `*`(I))), `*`(Physics:-Psigma[2]))), `-`(`*`(`/`(263697596812424, 996451875), `*`(Physics:-Psigma[1]))), `*`(`...
Typesetting:-mprintslash([C[10] := `+`(`*`(`*`(`/`(2185211616689851230363020476, 4204571658549609375), `*`(I)), `*`(Physics:-Psigma[2])), `*`(`/`(3226624781090887605597040906, 21022858292748046875), `... (222)
 

Notes: with the material above you can compute higher order values of C[n]. For that you need: 

  • Unassign C, as done above in two cases, to avoid interference of the results just computed.
 

  • Indicate more summation indices in the sequence summation_indices in summation_indices := n, m, k, l, p, q, r, s, t, u, as many as the maximum value of n in C[n].
 

  • Have in mind that the growth in size and complexity is significant, with each C[n] taking significantly more time than the computation of all the previous ones.
 

  • Re-execute the input line H := unapply(Formula(A, B, C, 10), N) and the loop for j to 9 do C[`+`(j, 1)] := solve(value(H(j)), C[`+`(j, 1)]) end do; 1.
 

Multivariable Taylor series of expressions involving anticommutative (Grassmannian) variables 

 The Physics:-Gtaylor command, for computing Taylor series of expressions involving anticommutative variables, got rewritten, now as a multivariable Taylor series command (same difference as in between the taylor and mtaylor commands) and combining two different approaches to handle, when possible, the presence of noncommutative and anticommutative variables. One is the standard approach for multi-variable expansions, it requires that the derivative of each function entering the expression being expanded commutes with the function itself. The second approach, for expansions with respect to anticommutative variables, separates the function into a "Body" and a "Soul", as is standard in supermathematics. 

Consider a set of anticommutative Grassmann variables theta__i with i = 1 .. n, forming a basis enlarged by their products, that satisfies `#mrow(mi(, so that The elements of the algebra involving these variables are linear combinations of the form 

 

 

 

where the coefficients `^`(`a__ij ...`, m)are complex numbers, `^`(a, 0) is called the body and everything else, `+`(a, `-`(`^`(a, 0))) = nil(a), is called the soul (nilpotent part). 

 

Consider now a mapping (function) f(a) which is assumed to be differentiable. Then f(a) can be defined by its Taylor expansion, 

 

 

 

That is the expansion computed by Gtaylor in Maple 2019. 

 

Examples 

 

> restart; 1; with(Physics); -1; Setup(anticommutativeprefix = {Theta, theta}, quantumoperators = {Theta, a})
 

[anticommutativeprefix = {Theta, theta}, quantumoperators = {Theta, a}] (223)
 

Consider 

> `/`(1, `*`(`+`(`*`(theta__1, `*`(theta__2)), theta__1, 1)))
 

Typesetting:-mprintslash([Physics:-`^`(`+`(1, Physics:-`*`(theta__1, theta__2), theta__1), -1)], [Physics:-`^`(`+`(1, Physics:-`*`(theta__1, theta__2), theta__1), -1)]) (224)
 

Step by step, `/`(1, `*`(`+`(`*`(theta__1, `*`(theta__2)), theta__1, 1))) is the application of the mapping  

> f := proc (u) options operator, arrow; `/`(1, `*`(u)) end proc
 

Typesetting:-mprintslash([f := proc (u) options operator, arrow; Physics:-`^`(u, -1) end proc], [proc (u) options operator, arrow; Physics:-`^`(u, -1) end proc]) (225)
 

to the element of the algebra  

> a := `+`(`*`(theta__1, `*`(theta__2)), theta__1, 1)
 

Typesetting:-mprintslash([a := `+`(1, Physics:-`*`(theta__1, theta__2), theta__1)], [`+`(1, Physics:-`*`(theta__1, theta__2), theta__1)]) (226)
 

The body of a is 

> body := eval(a, {theta__1 = 0, theta__2 = 0})
 

Typesetting:-mprintslash([body := 1], [1]) (227)
 

The soul of a is 

> soul := `+`(a, `-`(body))
 

Typesetting:-mprintslash([soul := `+`(Physics:-`*`(theta__1, theta__2), theta__1)], [`+`(Physics:-`*`(theta__1, theta__2), theta__1)]) (228)
 

and in view of 

> `*`(`^`(soul, 2)); -1; % = Expand(%)
 

Typesetting:-mprintslash([Physics:-`^`(`+`(Physics:-`*`(theta__1, theta__2), theta__1), 2) = 0], [Physics:-`^`(`+`(Physics:-`*`(theta__1, theta__2), theta__1), 2) = 0]) (229)
 

for the expression `/`(1, `*`(`+`(`*`(theta__1, `*`(theta__2)), theta__1, 1))), the Taylor expansion, mentioned in the introductory paragraph is 

> `+`(f(body), `*`((D(f))(body), `*`(soul)))
 

Typesetting:-mprintslash([`+`(1, `-`(Physics:-`*`(theta__1, theta__2)), `-`(theta__1))], [`+`(1, `-`(Physics:-`*`(theta__1, theta__2)), `-`(theta__1))]) (230)
 

The same computation, all in one go: 

> Gtaylor(Physics:-`^`(`+`(1, Physics:-`*`(theta__1, theta__2), theta__1), -1))
 

Typesetting:-mprintslash([`+`(1, `-`(Physics:-`*`(theta__1, theta__2)), `-`(theta__1))], [`+`(1, `-`(Physics:-`*`(theta__1, theta__2)), `-`(theta__1))]) (231)
 

  • Gauss integrals in this domain, that is, integrals over the exponential of a quadratic form of Grassmann variables, yield the determinant of the coefficient matrix of the quadratic form. These integrals play an important role in applications of anticommutative variables.
 

 

Problem. Taking into account that, with regards to Grassmann variables, differentiation and integration are the same operation, recover the determinant of the coefficient matrix with dimension N = 3 

 

Define the coefficient matrix and construct the exponential 

> N := 3; -1
 

> M := Matrix(N, symbol = m)
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (232)
 

> Theta__1 := Vector[row](N, proc (n) options operator, arrow; theta[`+`(n, 3)] end proc)
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (233)
 

> Theta__2 := Vector[column](N, proc (n) options operator, arrow; theta[n] end proc)
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-msub(Typesetting:-mi( (234)
 

> exp(Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(Theta__1, M), Theta__2))
 

Typesetting:-mprintslash([exp(`+`(`-`(Physics:-`*`(theta[1], `+`(`*`(m[1, 1], `*`(theta[4])), `*`(m[2, 1], `*`(theta[5])), `*`(m[3, 1], `*`(theta[6]))))), `-`(Physics:-`*`(theta[2], `+`(`*`(m[1, 2], `... (235)
 

The integral of this exponential can thus be obtained performing a multivariable Taylor series expansion, then differentiating (equivalent to integrating) with respect to the six theta[i] variables. 

> Theta := [theta[1], theta[2], theta[3], theta[4], theta[5], theta[6]]; -1
 

 

To avoid the default behavior of discarding terms of order 6 or higher in Theta, as in the other series commands of the Maple system, indicate the order term to be O(`*`(`^`(theta[i], 7))). 

> Gtaylor(exp(`+`(`-`(Physics:-`*`(theta[1], `+`(`*`(m[1, 1], `*`(theta[4])), `*`(m[2, 1], `*`(theta[5])), `*`(m[3, 1], `*`(theta[6]))))), `-`(Physics:-`*`(theta[2], `+`(`*`(m[1, 2], `*`(theta[4])), `*`...
 

Typesetting:-mprintslash([`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6...
Typesetting:-mprintslash([`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6...
Typesetting:-mprintslash([`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6...
Typesetting:-mprintslash([`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6...
Typesetting:-mprintslash([`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6...
(236)
 

Perform now the integration of the expanded exponential 

> diff(`+`(1, `-`(`*`(m[2, 2], `*`(Physics:-`*`(theta[2], theta[5])))), `-`(`*`(m[1, 2], `*`(Physics:-`*`(theta[2], theta[4])))), `-`(`*`(m[3, 1], `*`(Physics:-`*`(theta[1], theta[6])))), `-`(`*`(m[2, 1...
 

`+`(`*`(`+`(`*`(m[1, 1], `*`(m[2, 2])), `-`(`*`(m[1, 2], `*`(m[2, 1])))), `*`(m[3, 3])), `*`(m[3, 1], `*`(m[1, 2], `*`(m[2, 3]))), `*`(m[1, 3], `*`(`+`(`*`(m[2, 1], `*`(m[3, 2])), `-`(`*`(m[2, 2], `*`... (237)
 

Compare with the determinant of the coefficient matrix M 

> Delta := LinearAlgebra:-Determinant(M)
 

Typesetting:-mprintslash([Delta := `+`(`*`(m[1, 1], `*`(m[2, 2], `*`(m[3, 3]))), `-`(`*`(m[1, 1], `*`(m[2, 3], `*`(m[3, 2])))), `-`(`*`(m[1, 2], `*`(m[2, 1], `*`(m[3, 3])))), `*`(m[1, 2], `*`(m[2, 3],... (238)
 

> normal(`+`(Delta, `-`(`*`(`+`(`*`(m[1, 1], `*`(m[2, 2])), `-`(`*`(m[1, 2], `*`(m[2, 1])))), `*`(m[3, 3]))), `-`(`*`(m[3, 1], `*`(m[1, 2], `*`(m[2, 3])))), `-`(`*`(m[1, 3], `*`(`+`(`*`(m[2, 1], `*`(m[3...
 

0 (239)
 

  • An example with fermionic annihilation and creation operators, typical from quantum field theory.
 

 

> Setup(additionally, anticommutativeprefix = {Lambda, lambda}, op = Lambda)
 

 

 

`*`(`* Partial match of  '`, `*`(op, `*`(`' against keyword '`, `*`(quantumoperators, `*`(`' `)))))
_______________________________________________________
[anticommutativeprefix = {Lambda, Theta, lambda, theta}, quantumoperators = {Lambda, Theta, a}] (240)
 

 

Define corresponding annihilation and creation operators 

> am := Annihilation(Lambda)
 

Typesetting:-mprintslash([am := `#msup(mi( (241)
 

> ap := Creation(Lambda)
 

Typesetting:-mprintslash([ap := `#msup(mi( (242)
 

Note that the anticommutator of these two operators is equal to 1 (so they don't anticommute) 

> (%AntiCommutator = AntiCommutator)(am, ap)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (243)
 

while they do anticommute with all Grassmann variables 

> (%AntiCommutator = AntiCommutator)(am, theta)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (244)
 

> (%AntiCommutator = AntiCommutator)(ap, lambda)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (245)
 

and, because of Pauli's exclusion principle for fermions, the square of `#msup(mi( and `#msup(mi(are equal to 0.  

 

Consider now the product of exponentials 

> `*`(exp(`*`(theta, `*`(am))), `*`(exp(`+`(`-`(`*`(lambda, `*`(ap)))))))
 

Typesetting:-mprintslash([Physics:-`*`(exp(Physics:-`*`(theta, `#msup(mi( (246)
 

Because the corresponding exponents `*`(theta, `*`(`#msup(mi( and `+`(`-`(`*`(lambda, `*`(`#msup(mi( commute with their commutator, this product of exponentials can be combined using Hausdorff's formula 

> Physics:-`*`(exp(Physics:-`*`(theta, `#msup(mi(
 

Typesetting:-mprintslash([Physics:-`*`(exp(Physics:-`*`(theta, `#msup(mi( (247)
 

Both sides can be expanded in multivariable Taylor series, this time including the annihilation and creation operators as series variables since their commutation rules are all known ((%AntiCommutator = AntiCommutator)(am, ap), (%AntiCommutator = AntiCommutator)(am, theta), (%AntiCommutator = AntiCommutator)(ap, lambda)) and their square is equal to 0: 

> Gtaylor(Physics:-`*`(exp(Physics:-`*`(theta, `#msup(mi(
 

Typesetting:-mprintslash([`+`(1, `-`(Physics:-`*`(lambda, `#msup(mi( (248)
 

On the left-hand side is the product of the expansion of each exponential while on the right-hand side it is the expansion of the single exponential (no-trivially) combined taking into account the commutator of the exponents on the left-hand side. Verify that both expansions are one and the same: 

>
 

true (249)
 

  • Both Gtaylor and diff now handle Dagger(theta) as a differentiation or series variable in equal footing as theta itself
 

> exp(`+`(`*`(Dagger(theta), `*`(`#msup(mi(
 

Typesetting:-mrow(Typesetting:-mi( (250)
 

> exp(`+`(Physics:-`*`(Physics:-Dagger(theta), `#msup(mi(
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (251)
 

> diff(exp(`+`(Physics:-`*`(Physics:-Dagger(theta), `#msup(mi(
 

Typesetting:-mprintslash([`+`(`-`(`/`(1, 2)), Physics:-`*`(`#msup(mi( (252)
 

  • Depending on the case, multivariable Taylor expansions in the presence of not commutative variables can also be computed for functions of more than one variable and for unknown functions
 

> Gtaylor(F(x, theta, lambda), order = 3)
 

Typesetting:-mprintslash([`+`(F(0, 0, 0), `*`(x, `*`((D[1](F))(0, 0, 0))), `*`(`/`(1, 2), `*`(`^`(x, 2), `*`((D[1, 1](F))(0, 0, 0)))), `*`(`+`(`*`(x, `*`((D[1, 3](F))(0, 0, 0))), (D[3](F))(0, 0, 0)), ...
Typesetting:-mprintslash([`+`(F(0, 0, 0), `*`(x, `*`((D[1](F))(0, 0, 0))), `*`(`/`(1, 2), `*`(`^`(x, 2), `*`((D[1, 1](F))(0, 0, 0)))), `*`(`+`(`*`(x, `*`((D[1, 3](F))(0, 0, 0))), (D[3](F))(0, 0, 0)), ...
(253)
 

> Setup(noncommutativeprefix = {A, B})
 

[noncommutativeprefix = {A, B}] (254)
 

> Gtaylor(exp(`*`(`+`(A, B), `*`(delta))), delta = 0, order = 3)
 

`+`(1, `*`(delta, `*`(`+`(A, B))), `*`(`/`(1, 2), `*`(`^`(delta, 2), `*`(Physics:-`^`(`+`(A, B), 2))))) (255)
 

>
 

New SortProducts command 

A new command, SortProducts, receives an expression involving products and sorts the operands of these products according to the ordering indicated as the second argument, a list containing some or all of the operands of the product(s) found in the expression. The sorting of operands performed automatically takes into account any algebra rules set using Setup.  

 

Examples 

> restart; 1
 

> with(Physics); -1
 

Consider the product of the commutative a, b, c, d 

> P := `*`(a, `*`(b, `*`(c, `*`(d))))
 

Typesetting:-mprintslash([P := `*`(a, `*`(b, `*`(c, `*`(d))))], [`*`(a, `*`(b, `*`(c, `*`(d))))]) (256)
 

Reorder the operands c and b "in place". 

> SortProducts(P, [c, b])
 

`*`(a, `*`(c, `*`(b, `*`(d)))) (257)
 

Sort the operands b and c and put them to the left, then to the right 

> SortProducts(P, [c, b], totheleft)
 

`*`(c, `*`(b, `*`(a, `*`(d)))) (258)
 

> SortProducts(P, [c, b], totheright)
 

`*`(a, `*`(d, `*`(c, `*`(b)))) (259)
 

Set a prefix identifying noncommutative variables and related algebra rules for some of them, such that Z1, Z2, Z3 and Z4 commute between themselves, but none of them commute with Z5. 

> Setup(noncommutativeprefix = Z, %Commutator(Z1, Z2) = 0, %Commutator(Z1, Z3) = 0, %Commutator(Z2, Z3) = 0, %Commutator(Z3, Z4) = 0, %Commutator(Z2, Z4) = 0, %Commutator(Z1, Z4) = 0)
Setup(noncommutativeprefix = Z, %Commutator(Z1, Z2) = 0, %Commutator(Z1, Z3) = 0, %Commutator(Z2, Z3) = 0, %Commutator(Z3, Z4) = 0, %Commutator(Z2, Z4) = 0, %Commutator(Z1, Z4) = 0)
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (260)
 

> P := `*`(Z1, `*`(Z2, `*`(Z3, `*`(Z4))))
 

Typesetting:-mprintslash([P := Physics:-`*`(Z1, Z2, Z3, Z4)], [Physics:-`*`(Z1, Z2, Z3, Z4)]) (261)
 

Sort Z2 and Z3 "in place". 

> SortProducts(P, [Z3, Z2])
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z3, Z2, Z4)], [Physics:-`*`(Z1, Z3, Z2, Z4)]) (262)
 

> SortProducts(P, [Z3, Z2], totheright)
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z4, Z3, Z2)], [Physics:-`*`(Z1, Z4, Z3, Z2)]) (263)
 

The value of the commutator between Z1 and Z5 is not known to the system, so by default they are not sorted: 

> SortProducts(`*`(Z1, `*`(Z5)), [Z5, Z1])
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z5)], [Physics:-`*`(Z1, Z5)]) (264)
 

Force their sorting using their commutator or anticommutator 

> SortProducts(`*`(Z1, `*`(Z5)), [Z5, Z1], usecommutator)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (265)
 

> expand(`+`(Physics:-`*`(Z5, Z1), Physics:-Commutator(Z1, Z5)))
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z5)], [Physics:-`*`(Z1, Z5)]) (266)
 

> SortProducts(`*`(Z1, `*`(Z5)), [Z5, Z1], useanticommutator)
 

Typesetting:-mrow(Typesetting:-mo( (267)
 

> expand(`+`(`-`(Physics:-`*`(Z5, Z1)), Physics:-AntiCommutator(Z1, Z5)))
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z5)], [Physics:-`*`(Z1, Z5)]) (268)
 

Enter the product of P with Z5 at the end (so, to the right of P) and sort it with Z5 to the left. This is a case where Z5 does not commute with any of the other operands. Compare the results with and without the option usecommutator, 

> P := `*`(P, `*`(Z5))
 

Typesetting:-mprintslash([P := Physics:-`*`(Z1, Z2, Z3, Z4, Z5)], [Physics:-`*`(Z1, Z2, Z3, Z4, Z5)]) (269)
 

> SortProducts(P, [Z5, Z1])
 

Typesetting:-mprintslash([Physics:-`*`(Z1, Z2, Z3, Z4, Z5)], [Physics:-`*`(Z1, Z2, Z3, Z4, Z5)]) (270)
 

> SortProducts(P, [Z5, Z1], usecommutator)
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (271)
 

> SortProducts(P, [Z5, Z3], usecommutator)
 

Typesetting:-mrow(Typesetting:-mi( (272)
 

Documentation: "Physics updates", "A complete guide for performing tensor computations" and the "Mini-Course: Computer Algebra for Physicists" 

  • Prior to Maple 2019, the Physics package documentation contained

  • Help pages for each Physics command

  • A page on Physics conventions

  • Another page with Examples in different areas of physics

  • The "What's new in Physics" of each release with illustrations only shown there.

  • A number of Mapleprimes post describing the Physics project and showing how to use the package to tackle different problems.

  • Although this set is thorough in the information provided, the information is scattered. To address this situation, for Maple 2019, a single page, "Physics Updates" organizes and presents all those elusive links from b. to e. in one place, with hyperlinks to all the contents, and this page is now linked in all the Physics commands' help pages. Comments on practical ways to improve this presentation of information are welcome.

  • Likewise, one frequently asked question is on how to perform computations with tensors using the Physics packages, including the ones done in the past using GRTensor, a package developed in the 90's, not by Maplesoft, and which was at that time the standard for computer algebra and relativity. To address these question a new help page A Complete Guide for performing tensor computations using Physics, as an e-Handbook, is now part of the help system. This guide has three sections. Part I is all about tensors and their use in Euclidean spaces, Special Relativity, Quantum Mechanics and Classical field theory. Part II is all devoted to General Relativity. Part III is about transformations of coordinates on tensorial expressions.

  • The Physics package provides the mathematical tools for computations in several different areas. Maple users have mentioned other times that in a case like this the help pages themselves are not enough, that a tutorial approach would be an appropriate complement. In Maple 2019 we are thus adding a Mini-Course: Computer Algebra for Physicists. This course can be used as a tutorial, with 10 sections to be covered in 5 hands-on guided experiences of 2 hours each. The first part, 5 sections, is about Maple 101, for people who never used Computer Algebra, while the remaining 5 sections are all about using the Physics package.
 

Simplification of tensors, Pauli and Dirac matrices and KroneckerDelta 

 

  • Significant enhancement of the simplification of tensorial expressions
 

A new algorithm for normalizing tensorial expressions got implemented, following the ideas presented in the paper by L. R. U. Manssur, R. Portugal, and B. F. Svaiter,
Group-Theoretic Approach for Symbolic Tensor Manipulation, International Journal of Modern Physics C, Vol. 13, No. 07, pp. 859-879 (2002). 

 

Examples 

> restart; 1; with(Physics); -1; Setup(spacetimeindices = lowercaselatin)
 

[spacetimeindices = lowercaselatin] (273)
 

Define the following tensors, with no particular symmetry 

> Define(A, B, F, H, J, quiet)
 

{A, B, F, H, J, Physics:-Dgamma[a], Physics:-Psigma[a], Physics:-d_[a], Physics:-g_[a, b], Physics:-LeviCivita[a, b, c, d]} (274)
 

These tensors, A, B, F, H, and J respectively depend on 1 to 5 indices in the following tensorial expressions, all of which are actually equal to 0. The simplifier in Maple 2019 detects that fact by rewriting these expressions in the canonical / normal form explained in the reference mentioned above. To input the expressions with the contravariant indices as superscripts, type them preceded by ~, then right-click the input expression and select 2-D Math → Convert To → 2-D Math Input. 

>
 

> Simplify(e__1)
 

0 (275)
 

>
 

> Simplify(e__2)
 

0 (276)
 

>
 

> Simplify(e__3)
 

0 (277)
 

>
 

> Simplify(e__4)
 

0 (278)
 

The following example is less simple, define a tensor with three indices, that is symmetric with respect to the last two of them 

> Define(T[a, b, c], symmetric = [b, c])
 

 

`Defined objects with tensor properties`
{A[`~c`], B[a, `~d`], Physics:-Dgamma[a], F[e, d, `~j`], H[a, f, `~f`, j], J[e, c, j, k, f], Physics:-Psigma[a], T[a, b, c], Physics:-d_[a], Physics:-g_[a, b], Physics:-LeviCivita[a, b, c, d]} (279)
 

So 

> `+`(T[1, 2, 3], `-`(T[1, 3, 2]))
 

`+`(T[1, 2, 3], `-`(T[1, 3, 2])) (280)
 

> Simplify(`+`(T[1, 2, 3], `-`(T[1, 3, 2])))
 

0 (281)
 

If we now swap a with c and take the difference we get an antisymmetric tensorial expression 

> TT := `+`(`*`(2, `*`(Antisymmetrize(T[a, b, c], [a, c]))))
 

Typesetting:-mprintslash([TT := `+`(T[a, b, c], `-`(T[c, b, a]))], [`+`(T[a, b, c], `-`(T[c, b, a]))]) (282)
 

So, by construction, the following is equal to 0 even when none of the terms is; detecting situations like this one is part of the intrinsic efficiency of the group theoretic approach 

> expand(`*`(A[a], `*`(A[c], `*`(TT)))) = 0
 

`+`(`*`(A[`~a`], `*`(A[`~c`], `*`(T[a, b, c]))), `-`(`*`(A[`~a`], `*`(A[`~c`], `*`(T[c, b, a]))))) = 0 (283)
 

> Simplify(`+`(`*`(A[`~a`], `*`(A[`~c`], `*`(T[a, b, c]))), `-`(`*`(A[`~a`], `*`(A[`~c`], `*`(T[c, b, a]))))) = 0)
 

0 = 0 (284)
 

  • The Pauli and Dirac matrices are implemented as 4-vectors
 

 

In Maple 2019, to compute the matrix components use Library:-RewriteInMatrixForm and Library:-PerformMatrixOperations 

> restart; 1; with(Physics); -1; Setup(spaceindices = lowercaselatin)
 

[spaceindices = lowercaselatin] (285)
 

Since Pauli matrices are now defined as a 4-vector, all the keywords for tensors automatically work 

> Psigma[definition]
 

Typesetting:-mprintslash([Typesetting:-mcomplete(Typesetting:-msub(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (286)
 

Likewise, you can visualize tensor components the usual way using TensorArray 

> TensorArray([%Commutator(Physics:-Psigma[a], Physics:-Psigma[b]) = `*`(`*`(2, `*`(I)), `*`(Physics:-LeviCivita[a, b, `~c`], `*`(Physics:-Psigma[c]))), %AntiCommutator(Physics:-Psigma[a], Physics:-Psig...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (287)
 

To see the matrix expression of these commutators and anticommutators of Pauli matrices use the option performmatrixoperations. For example, for the first block of identities involving commutators, 

> TensorArray([%Commutator(Physics:-Psigma[a], Physics:-Psigma[b]) = `*`(`*`(2, `*`(I)), `*`(Physics:-LeviCivita[a, b, `~c`], `*`(Physics:-Psigma[c]))), %AntiCommutator(Physics:-Psigma[a], Physics:-Psig...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (288)
 

The simplifier now knows more about Pauli matrices 

> `+`(`*`(2, `*`(Psigma[1], `*`(Psigma[2]))))
 

Typesetting:-mprintslash([`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[2]), Physics:-`*`(Physics:-Psigma[2], Physics:-Psigma[1]))], [`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[2]), P... (289)
 

> Simplify(`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[2]), Physics:-`*`(Physics:-Psigma[2], Physics:-Psigma[1])))
 

0 (290)
 

> `*`(`^`(Psigma[2], 2))
 

Physics:-`^`(Physics:-Psigma[2], 2) (291)
 

> Simplify(Physics:-`^`(Physics:-Psigma[2], 2))
 

1 (292)
 

> `*`(Psigma[1], `*`(Psigma[3]))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[3])], [Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[3])]) (293)
 

> Simplify(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[3]))
 

`+`(`-`(`*`(`+`(I), `*`(Physics:-Psigma[2])))) (294)
 

These two library routines are the ones used to rewrite tensorial expressions in matrix form or to perform the corresponding matrix operations 

> Library:-RewriteInMatrixForm(`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[2]), Physics:-`*`(Physics:-Psigma[2], Physics:-Psigma[1])))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (295)
 

> Library:-PerformMatrixOperations(`+`(Physics:-`*`(Physics:-Psigma[1], Physics:-Psigma[2]), Physics:-`*`(Physics:-Psigma[2], Physics:-Psigma[1])))
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (296)
 

The same works for the Dirac matrices 

> Dgamma[definition]
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (297)
 

> TensorArray(%AntiCommutator(Physics:-Dgamma[`~mu`], Physics:-Dgamma[`~nu`]) = `+`(`*`(2, `*`(Physics:-g_[`~mu`, `~nu`]))))
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (298)
 

New in Maple 2019, when Physics is loaded the standard representation of Dirac matrices is automatically loaded too, corresponding to the contravariant components of gamma[`~mu`] 

> Dgamma[`~`]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (299)
 

>
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (300)
 

The definition of gamma[5] is also visible using the keyword definition 

> Dgamma[5, definition]
 

Typesetting:-mprintslash([Physics:-Dgamma[5] = Physics:-Dgamma[`~5`], Physics:-Dgamma[`~5`] = `+`(`-`(`*`(`+`(I), `*`(Physics:-`*`(Physics:-Dgamma[`~0`], Physics:-Dgamma[`~1`], Physics:-Dgamma[`~2`], ... (301)
 

Verify the first three of these identities 

> value(TensorArray([Physics:-Dgamma[5] = Physics:-Dgamma[`~5`], Physics:-Dgamma[`~5`] = `+`(`-`(`*`(`+`(I), `*`(`*`(Physics:-Dgamma[`~0`], Physics:-Dgamma[`~1`], Physics:-Dgamma[`~2`], Physics:-Dgamma[...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn( (302)
 

For the fourth identity 

> (Physics:-Dgamma[5] = Physics:-Dgamma[`~5`], Physics:-Dgamma[`~5`] = `+`(`-`(`*`(`+`(I), `*`(`*`(Physics:-Dgamma[`~0`], Physics:-Dgamma[`~1`], Physics:-Dgamma[`~2`], Physics:-Dgamma[`~3`]))))), `*`(Ph...
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (303)
 

> value(TensorArray([Physics:-Dgamma[5] = Physics:-Dgamma[`~5`], Physics:-Dgamma[`~5`] = `+`(`-`(`*`(`+`(I), `*`(`*`(Physics:-Dgamma[`~0`], Physics:-Dgamma[`~1`], Physics:-Dgamma[`~2`], Physics:-Dgamma[...
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mfenced(Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetti... (304)
 

You can compute with the tensor components and later represent them in matrix form, or perform the corresponding matrix operations 

> `+`(`*`(Dgamma[1], `*`(Dgamma[2])), Dgamma[0])
 

Typesetting:-mprintslash([`+`(Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[2]), Physics:-Dgamma[4])], [`+`(Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[2]), Physics:-Dgamma[4])]) (305)
 

> Library:-RewriteInMatrixForm(`+`(Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[2]), Physics:-Dgamma[4]))
 

Typesetting:-mrow(Typesetting:-mrow(Typesetting:-mi( (306)
 

> Library:-PerformMatrixOperations(`+`(Physics:-`*`(Physics:-Dgamma[1], Physics:-Dgamma[2]), Physics:-Dgamma[4]))
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mrow(Typesetting:-mn( (307)
 

 

> `*`(`^`(Dgamma[mu], 2), `*`(Dgamma[`~lambda`], `*`(Dgamma[`~nu`], `*`(Dgamma[`~rho`], `*`(Dgamma[`~sigma`])))))
 

Typesetting:-mprintslash([Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~mu`])], [Physics:-... (308)
 

> Simplify(Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~mu`]))
 

Typesetting:-mprintslash([`+`(`*`(2, `*`(Physics:-`*`(Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`]))), `*`(2, `*`(Physics:-`*`(Physics:-Dgamma... (309)
 

> TensorArray(`+`(Physics:-`*`(Physics:-Dgamma[mu], Physics:-Dgamma[`~lambda`], Physics:-Dgamma[`~nu`], Physics:-Dgamma[`~rho`], Physics:-Dgamma[`~sigma`], Physics:-Dgamma[`~mu`]), `-`(`*`(2, `*`(Physic...
 

Typesetting:-msub(Typesetting:-mi( (310)
 

>
 

Typesetting:-msub(Typesetting:-mi( (311)
 

  • In Maple 2019, KroneckerDelta is not a tensor
 

 

Although in an Euclidean space the Kronecker delta symbol is a tensor, its components do not change under a transformation of coordinates, that is not the case in a Minkowski or curved spacetime. Also, KroneckerDelta is more often than otherwise used taking delta[a, b] = 1 when a = b, while due to Einstein's sum rule for repeated indices, if KroneckerDelta were a tensor and a is a tensor index, then delta[a, a] = `*`(dimension, `*`(of, `*`(space))). To avoid this ambiguity of notation, in Maple 2019 KroneckerDelta is not implemented as a tensor, but as the standard non-tensorial Kronecker delta symbol.

For the cases where you need to use it as a tensor, for example when entering commutation relations in quantum mechanics using tensorial notation, you can either use the metric itself g[a, b], or define a Kronecker delta tensor for that purpose. For example
 

> restart; 1; with(Physics); -1; Setup(spacetimeindices = lowercaselatin, metric = Euclidean, dimension = 3)
 

 

 

 

`The dimension and signature of the tensor space are set to `[3, `- - +`]
`The Euclidean metric in cartesian coordinates`
`*`(`Changing the signature of the tensor spacetime to: `, `*`(`+ + +`))
[dimension = 3, metric = {(1, 1) = 1, (2, 2) = 1, (3, 3) = 1}, spacetimeindices = lowercaselatin] (312)
 

This does not return the dimension of space 

> KroneckerDelta[a, a]
 

1 (313)
 

This returns the dimension of space: 

> g_[a, a]
 

3 (314)
 

You can define a Kronecker delta; tensor using Define with ease, for example defining the covariant components of the tensor as follows 

> delta[j, k] = Matrix(3, proc (j, k) options operator, arrow; KroneckerDelta[j, k] end proc)
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (315)
 

>
 

 

`Defined objects with tensor properties`
{Physics:-Dgamma[a], Physics:-Psigma[a], Physics:-d_[a], delta[j, k], Physics:-g_[a, b], Physics:-LeviCivita[a, b, c]} (316)
 

Now you have 

> delta[]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (317)
 

> delta[a, `~b`, matrix]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (318)
 

> delta[`~`]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (319)
 

And the trace: 

> delta[trace]
 

3 (320)
 

So this is not equal to 1 

> delta[a, a]
 

delta[a, a] (321)
 

> SumOverRepeatedIndices(delta[a, a])
 

3 (322)
 

This delta tensor, well defined in an Euclidean space, however changes when the space is not Euclidean. For example: 

> Setup(signature = `-`)
 

[signature = `- - +`] (323)
 

> delta[]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (324)
 

> delta[a, `~b`, matrix]
 

Typesetting:-mrow(Typesetting:-msub(Typesetting:-mi( (325)
 

> delta[trace]
 

-1 (326)