-
Notifications
You must be signed in to change notification settings - Fork 0
/
nonrelativistic_qm_ops.mac
209 lines (167 loc) · 6.21 KB
/
nonrelativistic_qm_ops.mac
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
/* Maxima code for working with quantum mechanical operators
Barton Willis
Professor of Mathematics
University of Nebraska at Kearney
Copyright (c) Barton Willis, 2022
GNU GENERAL PUBLIC LICENSE Version 3, 29 June 2007
*/
/* Make sure that "." is noncommutative */
remove(".",commutative);
commutes_with[a,b]:= false; /*backstop for undefined */
/* Declare the reduced Planck constant to be constant. That makes ħ outative for a linear operator.*/
declare(ħ,constant);
assume(ħ > 0);
/* Declare position and momentum operators. */
declare([Qx,Qy,Qz,Px,Py,Pz],linear);
declare([Qx,Qy,Qz,Px,Py,Pz],operator);
/* Predicates for detecting these operators */
Qx_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Qx);
Qy_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Qy);
Qz_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Qz);
Px_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Px);
Py_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Py);
Pz_p(e) := block([inpart : true], not mapatom(e) and inpart(e,0) = 'Pz);
/* Give them formulas and adjoints*/
put(Px, lambda([q], -%i*ħ*diff(q,x)), 'formula);
put(Px, Px, 'operator_adjoint);
put(Py, lambda([q], -%i*ħ*diff(q,y)), 'formula);
put(Py, Py, 'operator_adjoint);
put(Pz, lambda([q], -%i*ħ*diff(q,z)), 'formula);
put(Pz, Pz, 'operator_adjoint);
put(Qx, lambda([q], x*q), 'formula);
put(Qx, Qx, 'operator_adjoint);
put(Qy, lambda([q], y*q), 'formula);
put(Qy, Qy, 'operator_adjoint);
put(Qz, lambda([q], z*q), 'formula);
put(Qz, Qz, 'operator_adjoint);
/* Implement the rules
Qx Qy --> Qy Qx
Qx Qz --> Qz Qy
Qx σx --> σx Qx
Qx σy --> σy Qx
Qx σz --> σz Qx */
commutes_with[Qx,Qy] : true;
commutes_with[Qx,Qz] : true;
commutes_with[Qx,σx] : true;
commutes_with[Qx,σy] : true;
commutes_with[Qx,σz] : true;
simp_Qx(e) := block([inflag : true],
if mapatom(e) then simpfuncall('Qx,e)
elseif commutes_with[Qx,inpart(e,0)] then (
inpart(e,0)(apply('Qx, args(e))))
else simpfuncall(Qx,e));
simplifying('Qx, 'simp_Qx);
/* Implement the rules
Qy Qz --> Qz Qy
Qy σx --> σx Qy
Qy σy --> σy Qy
Qy σz --> σz Qy */
commutes_with[Qy,Qz] : true;
commutes_with[Qy,σx] : true;
commutes_with[Qy,σy] : true;
commutes_with[Qy,σz] : true;
simp_Qy(e) := block([inflag : true],
if mapatom(e) then simpfuncall('Qy,e)
elseif commutes_with[Qy,inpart(e,0)] then (
inpart(e,0)(apply('Qy, args(e))))
else simpfuncall(Qy,e));
simplifying('Qy, 'simp_Qy);
/* Implement the rules
Qz σx --> σx Qz
Qz σy --> σy Qz
Qz σz --> σz Qz */
commutes_with[Qz,σx] : true;
commutes_with[Qz,σy] : true;
commutes_with[Qz,σz] : true;
simp_Qz(e) := block([inflag : true],
if mapatom(e) then simpfuncall('Qz,e)
elseif commutes_with[Qz,inpart(e,0)] then (
inpart(e,0)(apply('Qz, args(e))))
else simpfuncall(Qz,e));
simplifying('Qz, 'simp_Qz);
/* Implement the rules
Px Qx --> Qx Px - %i*ħ
Px Qy --> Qy Px
Px Qz --> Qz Px
Px Py --> Py Px
Px Pz --> Pz Px
Px σx --> σx Px
Px σy --> σy Px
Px σz --> σz Px */
commutes_with[Px,Py] : true;
commutes_with[Px,Pz] : true;
commutes_with[Px,Qy] : true;
commutes_with[Px,Qz] : true;
commutes_with[Px,σx] : true;
commutes_with[Px,σy] : true;
commutes_with[Px,σz] : true;
simp_Px(e) := block([inflag : true],
if mapatom(e) then simpfuncall('Px,e)
elseif Qx_p(e) then Qx(Px(first(e))) - %i*ħ*first(e)
elseif commutes_with[Px,inpart(e,0)] then (
inpart(e,0)(apply('Px, args(e))))
else simpfuncall('Px,e));
simplifying('Px, 'simp_Px);
/*Pauli matrices */
declare([σx,σy,σz],linear, [σx,σy,σz], operator);
put(σx, lambda([q], σx.q), 'formula);
put(σy, lambda([q], σy.q), 'formula);
put(σz, lambda([q], σz.q), 'formula);
put(σx,σx,'operator_adjoint);
put(σy,σy,'operator_adjoint);
put(σz,σz,'operator_adjoint);
σx_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'σx);
σy_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'σy);
σz_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'σz);
put(σx, lambda([q], matrix([0,1],[1,0]).q), formula);
put(σy, lambda([q], matrix([0,-%i],[%i,0]).q), formula);
put(σz, lambda([q], matrix([1,0],[0,-1]).q), formula);
/* angular momentum */
/*Tell Maxima that Lx, Ly, and Lz are linear operators.*/
declare([Lx, Ly, Lz],linear, [Lx, Ly, Lz],operator);
Lx_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'Lx);
Ly_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'Ly);
Lz_p(e) := block([inflag : true], not mapatom(e) and inpart(e,0) = 'Lz);
/* Give Lx, Ly, and Lz adjoints. */
put(Lx, Lx, 'operator_adjoint);
put(Ly, Ly, 'operator_adjoint);
put(Lz, Lz, 'operator_adjoint);
/*Give Lx, Ly, and Lz formula */
put(Lx, lambda([q], -%i*ħ*(y * diff(q,z) - z*diff(q,y))), 'formula);
put(Ly, lambda([q], -%i*ħ*(z * diff(q,x) - x*diff(q,z))), 'formula);
put(Lz, lambda([q], -%i*ħ*(x * diff(q,y) - y*diff(q,x))), 'formula);
/* Implement the rules Lz Ly --> Ly Lz -%i ħ Lx and Lz Lx --> Lx Lz + %i ħ Ly */
simp_Lz(e) :=
if Ly_p(e) then Ly(Lz(first(e))) - %i*ħ*Lx(first(e))
elseif Lx_p(e) then Lx(Lz(first(e))) + %i*ħ*Ly(first(e))
else simpfuncall('Lz,e);
simplifying('Lz, 'simp_Lz);
/* Implement the rule Ly Lx --> Lx Ly - %i ħ Lz */
simp_Ly(e) :=
if Lx_p(e) then Lx(Ly(first(e))) - %i*ħ*Lz(first(e))
else simpfuncall('Ly,e);
simplifying('Ly, 'simp_Ly);
/* Implement the rules
σx^2 --> I,
σx σy --> σy σx - 2 %i σz,
σx σz --> σz σx - 2 %i σy. */
simp_σx(e) := block([inflag : true],
if σx_p(e) then first(e)
elseif σy_p(e) then σy(σx(first(e))) - 2*%i*σz(first(e))
elseif σz_p(e) then σz(σx(first(e))) - 2*%i*σy(first(e))
else simpfuncall('σx,e));
simplifying('σx, 'simp_σx);
/* Implement the rules
σy^2 --> I,
σy σz --> σz σy - 2*%i σx. */
simp_σy(e) := block([inflag : true],
if σy_p(e) then first(e)
elseif σz_p(e) then σz(σy(first(e))) - 2*%i*σx(first(e))
else simpfuncall('σy,e));
simplifying('σy, 'simp_σy);
/* Implement the rule
σz σz --> 1 */
simp_σz(e) := block([inflag : true],
if σz_p(e) then first(e)
else simpfuncall('σz,e));
simplifying('σz, 'simp_σz);