Cuantica: A Package for Quatum Mechanical Calculus

Lic. Renan Cabrera Lafuente
Academia Nacional de Ciencias de Bolivia

renanbo@hotmail.com
rencabla@ceibo.entelnet.bo
La Paz Bolivia -2000

Abstract
This is a package made for quantum mechanical calculus and  specially for symbolic purpose. We know that this field is very very broad, but I put in this firsts versions several important topics, that I know will be useful saving time from tedious calculus.

-Warning: This notebook requires the package "cuantica.m" version 1.2 and the package "versora.m" version 1.1 or above!!!

- "cuantica.m" is a package published at:
http://www.mathsource.com/Content/Applications/Physics/0211-374

Needs["cuantica`", 
      "c:\\renan\\mathematica\\book\\code\\cuantica.m"]

Needs["versora`","c:\\Renan\\mathematica\\Book\\code\\versora.m"]\(\(SetCoordinates[e]\)[\ 
Cartesian[{x, y, z}], \ 
        Cylindrical[{rhc, thc}]\ , \ Spherical[{r, th, ph}]\ \ \ ]\)

\(Format[ph] = φ; 
Format[th] = θ; \)

The Ket Hilbert Space

Introduction

First of all we have to know how to define Orthonormal Basis in a Hilbert space, This task is done by SetDiscreteOB which first argument is the Label at which is attached the particular Hilbert Space.

SetDiscreteOB[ Ket[Sz, {-1, 1, 2}]]\(SetDiscreteOB[\ Ket[a, {{2, 5}, {0, 3}}]]
\)

The inner product is represented by the HIP function, that operates between the Bra  and the Ket spaces as usual in Dirac notation

\(\(
HIP[\ Bra[Sz, 1/2], Ket[Sz, \(-1\)/2]\ \ ]\)\)

\(\(
HIP[p\ Bra[a, {3, 0}], \ 2\ Ket[a, {3, 0}]\ ]\)\)

HIP[ p ( Bra[a, b] +  Bra[1, 2]) , 
      2  Ket[c, d] +  λ Ket[u, v] ]

HO[A]HO[ 2  I Ket[a, b] + Ket[p, q]]

Transformation of Basis

Let us make some basic Quantum Mechanical Calculus in a 3 dimensional space. First lets define two Ortonornal Base.

BaseA  =  SetDiscreteOB[ Ket[A, {1, 3}]]BaseB  =  SetDiscreteOB[ Ket[B, {1, 3}]]

The function HEP gives the external product

HEP[Ket[1, 1], Bra[2, 3]]

In this way we can define the Identity Operator

\(IDA\  = \ 
      Evaluate[
		Inner[HEP, BaseA\ , HO[BaseA], Plus]

          	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]&
		\)\(IDB\  = \ 
      Evaluate[
		Inner[HEP, BaseB\ , HO[BaseB], Plus]

          	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]&\)

Now we can define a posible associated Hermitic Operator to the A's Basis

\(AOperator\  = \ 
      Evaluate[
	\ 
          HEP[
			\ \ A1\ Ket[A, 1], Bra[A, 1]] + \ 
            A2\ HEP[Ket[A, 2], Bra[A, 2]] + A3\ HEP[Ket[A, 3], Bra[A, 3]]\ 

          		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]&
\)\(Table[
HIP[\ \ Bra[A, j]\ , AOperator[Ket[A, i]]\ ]

      \ \ \ \ \ \ \ \ \ , {i, 1, 3}, {j, 1, 3}]\)

Next we define a Unitary Operator for transforming basis.

\(BaseB\  = \ SetDiscreteOB[\ Ket[B, {1, 3}]]; 


    UA2B\  = \ \ 
      Evaluate[
		Inner[HEP, BaseB\ , HO[BaseA], Plus]

          	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]&\)

And it is possible to verify

UA2B[ Ket[A, 1]]

Lets see its matrix form in the "A" Basis

\(UA2BMatrix\  = \ \ 
      Table[
HIP[\ \ Bra[A, i]\ , UA2B[Ket[A, j]]\ \ ]

        \ \ \ \ \ \ \ \ \ , {i, 1, 3}, {j, 1, 3}]\)

The Hermitic operator over the Unitary Operator gives its inverse as we expect

UB2A  =  HO[ UA2B ]

\(Evaluate[
UB2A[UA2B[#]]
]&\)

Let | ψ , 1 > be a particular ket in the 3D Hilber space. We can represent this ket
in the both the A's system and B's system.

ψA  =  IDA[Ket[ψ, 1]]\(ψB\  = \ IDB[Ket[ψ, 1]]
\)

Now the problem is how to relate their components in both Basis.

If we find the Matrix form of both representations

Table[  HIP[Bra[A, i], ψA]  ,  {i, 1, 3} ]  // TableForm

Table[  HIP[Bra[A, i], ψB]  ,  {i, 1, 3} ]  // TableForm

The later could be also found by multiplying the B's  representation by the Unitary matrix A - > B  in the A's represntatión.

UA2BMatrix. Table[  HIP[Bra[B, i], ψB]  ,  {i, 1, 3} ]  // 
      TableForm

Position Representation. Schrodinger Equation

Potential Well

Let V[ parameters ][x] the potential well

A particle in a potential well can be expressed by the following parameters:
L = Length of the well
V1 = Potential value. Region 1.
V3 = Potential value. Region 3.
ϵ  =  Energy of the particle.
m = mass of the particle.

V[L_, V1_, V3_][x_] :=  V1  /;  x < 0V[L_, V1_, V3_][x_] :=  0  /;  x < L  &&  x > 0\(\(V[L_, V1_, V3_]\)[x_]\  := \ V3\  /; \ x > L
\)

A particular case:

Plot[ V[6, 2, 3][x] , {x, -2, 8}, Frame -> True, 
      Axes -> False, PlotLabel -> "Potential Well", 
      FrameLabel -> {"x", "ϵ"}, RotateLabel -> False]

Case 1: Energy particle E < V1 and E< V3

The proposed solution is:

\(ψ1 = \ A\ Exp[\ α\ x]; 
	

    ψ2 = \ B\ Sin[k\ x\  + δ]; 
		

    ψ3 = \ C\ Exp[\(-β\)\ x]; 


    Toϵ\  = 
      \ {\ α\  -> \ \@\(\(2  m \((V1 - ϵ)\)\)\/ℏ\^2\), 
        k -> \@\(\(2  m\ ϵ\)\/ℏ\^2\), 
        β -> \ \@\(\(2  m \((V3 - ϵ)\)\)\/ℏ\^2\)}; \)

The parameters to fix are A, B,  C, δ and the energy ϵ . There are two conditions in each point of potential change (Continuity and Derivability). These leave us one parameter free wich will be used to normalize the wave function to 1. 
Instead of the continuity restriction it is better to take the Derivative of the Logaritm.

Now appliying the Derivative of the Logaritm and expressing α and β by k 
(all of them are functions of the energy ϵ).

\(\(\ 

    	{D[Log[ψ1], x] == D[Log[ψ2], x] /. \ x -> 0, 
	
        D[Log[ψ2], x] == D[Log[ψ3], x]	 /. x -> L} /. 
      \ {α\  -> \ \@\(\(2  m\ V1\)\/ℏ\^2 - k\^2\)\ \ , 
        β\  -> \ \@\(\(2  m\ V3\)\/ℏ\^2 - k\^2\)}
\)\)

And we get a system of two equations with two unknown variables  k and  δ.
Expressing the tangent function by the sine function

And regarding the periodiod π   of the cotangent function we obtain:

             \(TraditionalForm
      \`k\ L = π\ n - \(sin\^\(-1\)\)(\(ℏ\ k\)\/\@\(2\ m\ V1\)) - 
          \(sin\^\(-1\)\)(\(ℏ\ k\)\/\@\(2\ m\ V2\))\)

Parameters  = 
       {A -> 1, ℏ  ->  1 ,  m  -> 1 ,  V1  ->  2 ,  
        V3  -> 3 ,  L  ->  6};

And the following graphic shows the particular and discrete values of  k into the intersections. 
( Case V1 < V2 )

\(\(
Plot[
	
      Evaluate[{\ k\ L, \ 
            Table[n\ π - ArcSin[\(ℏ\ k\)\/\@\(2\ m\ V1\)] - 
                ArcSin[\(ℏ\ k\)\/\@\(2\ m\ V3\)], {n, 1, 5}]} /. 
          Parameters\ ]
, {k, 0, \@\(2\ m\ V1\)\/ℏ /. Parameters}, 
      Frame -> True, Axes -> False]\)\)

And we can guess the solutions

Solk  =  {0.45, 0.9, 1.35, 1.78};

Solϵ =  ((ℏ^2)/(2 m)) (Solk^2)  /. Parameters;

Now lets find the solution of  δ, B and C.

Solδ  =  
      ArcSin[((ℏ #)/(Sqrt[2 m V1]))] /. Parameters  &/@ Solk

In order to continue we have to define some parameters. We can chose A arbitrarily and leave the normalization changed temporary.

\(SolB\  = 
      Function[{δ}, \ 
          Evaluate[
            B /. \(Solve[
                  \(ψ1 == ψ2\  /. \ x -> 0\) /. Parameters\ , {B}]
                  \)[\([1]\)]\ ]

          		\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]/@\ Solδ\)

SolC  = 
      C /.  
                Solve[ψ2 == ψ3  /. { x -> L }, {C}][[1]]
                   /.  β  ->  
                (Sqrt[((2  m V3)/(ℏ^2)) - (k^2)]) /. Parameters /.  
          Map[ MapThread[Rule, {{k, δ, B}, #}]& , 
            Transpose[{Solk, Solδ, SolB}]] // Flatten

\(SolRule = 
      Transpose[{\(Outer[Rule, {ϵ}, Solϵ]\)[\([1]\)], 

          \(Outer[Rule, {δ}, Solδ]\)[\([1]\)], 

          \(Outer[Rule, {B}, SolB]\)[\([1]\)], 

          \(Outer[Rule, {C}, SolC]\)[\([1]\)]}]\)

ψ1List =  
      ψ1 /. Toϵ /. Parameters  /. SolRule ;

ψ2List  = 
      ψ2 /. Toϵ  /. SolRule /. Parameters;

ψ3List  = 
      ψ3 /. Toϵ /. SolRule /. Parameters;

\(\(
Show[

      Plot[\ Evaluate[ψ1List], {x, \(-2\), 0}, 
        DisplayFunction -> Identity]
, 

      Plot[\ Evaluate[ψ2List], {x, 0, L /. Parameters}, 
        DisplayFunction -> Identity]
, 

      Plot[\ Evaluate[ψ3List], {x, L /. Parameters, 
          L + 2 /. Parameters}, DisplayFunction -> Identity]
, 
	
      Plot[
	Evaluate[\ \(V[L, V1, V3]\)[x] /. Parameters]\ 

        		, {x, \(-1\), 7}, PlotStyle -> {Thickness[0.007]}, 
        DisplayFunction -> Identity]
, 
	
      Plot[\ Evaluate[Solϵ], {t, 0, L /. Parameters}, Frame -> True, 
        PlotStyle -> {Dashing[{0.01, 0.01}]}, DisplayFunction -> Identity]

      
, Axes -> False, Frame -> True, DisplayFunction -> $DisplayFunction, 
      PlotLabel -> "Quantum Potential Well", 
	
      FrameLabel -> {"x", "ϵ"}, RotateLabel -> False]



    \)\)

Case 2: Energy particle E > V1 and E< V2

We make the variable change ϵ ->  ((2 m ϵ)/(ℏ^2))

\(ψ1 = \ A1\ Exp[\ I\ k1\ x] + B1\ Exp[\(-I\)\ k1\ x]; 
	

    ψ2 = \ A2\ Exp[\ I\ k2\ x] + B2\ Exp[\(-I\)\ k2\ x]; 
		

    ψ3 = \ C\ Exp[\(-β\)\ x]; 


    kToϵ\  = 
      \ {\ k1\  -> \ \@\(ϵ - V1\), k2 -> \@ϵ, 
        β -> \ \@\(V3 - ϵ\)}; \)

Parameters  = 
       {C -> 1000, L ->  6 ,  V1  ->  2,  V3  ->  4, 
        ϵ -> 3};

\(\(	\({\(ψ1\  == \ ψ2\  /. \ x -> 0\) /. Parameters
	, 
	
      \(D[ψ1, x] == D[\ ψ2, x] /. \ x -> 0\) /. Parameters
	, 
	
      \(ψ2\  == \ ψ3\  /. \ x -> L\) /. Parameters
	, 
	
      \(D[\ ψ2, x] == D[ψ3, x]	 /. x -> L\)	 /. Parameters
	}; 
    \)\)\)

\(\(Sol = 
      Solve[	
	{ψ1\  == \ ψ2\  /. \ x -> 0
	, 
	
            D[ψ1, x] == D[\ ψ2, x] /. \ x -> 0
	, 
	
            ψ2\  == \ ψ3\  /. \ x -> L
	, 
	
            D[\ ψ2, x] == D[ψ3, x]	 /. x -> L	
	}, 

          	{A1, B1, A2, B2}] // Simplify; \)\)

\(ψ1Sol = \ 
      \(\(A1\ Exp[\ I\ k1\ x] + B1\ Exp[\(-I\)\ k1\ x] /. \ Sol\) /. 
          kToϵ\) /. Parameters; 
	

    ψ2Sol = \ 
      \(\(A2\ Exp[\ I\ k2\ x] + B2\ Exp[\(-I\)\ k2\ x] /. Sol\) /. 
          kToϵ\) /. Parameters; 
		

    ψ3Sol = \ 
      \(C\ Exp[\(-β\)\ x] /. kToϵ\) /. Parameters; \)

\(Show[
	
      Plot[\ Re[ψ1Sol]\ , {x, \(-12\), 0}, 
        DisplayFunction -> Identity\ ], 
		
      Plot[\ Re[ψ2Sol]\ , {x, 0, L /. Parameters}\ , 
        DisplayFunction -> Identity], 
		
      Plot[\ Re[ψ3Sol]\ , {x, L /. Parameters, 2  L /. Parameters}\ , 
        DisplayFunction -> Identity]
, 
	
      Plot[
	Evaluate[\ \(V[L, V1, V3]\)[x] /. Parameters]\ 

        		, {x, \(-12\), 12}, PlotStyle -> {Thickness[0.007]}, 
        DisplayFunction -> Identity]
, 
	
      Plot[\ ϵ\  /. Parameters, {t, \(-12\), L /. Parameters}, 
        Frame -> True, PlotStyle -> {Dashing[{0.01, 0.01}]}, 
        DisplayFunction -> Identity]
, Frame -> True, 
      DisplayFunction -> $DisplayFunction, PlotLabel -> "Re[ψ]", 
      FrameLabel -> {"x", "ϵ"}, RotateLabel -> False]\)

Potential Barrier

Clear[V]

V[L_, V2_][x_] :=  0  /;  x < 0V[L_, V2_][x_] :=  V2  /;  x < L  &&  x > 0\(\(V[L_, V2_]\)[x_]\  := \ 0\  /; \ x > L
\)

Plot[ V[6, 3][x] , {x, -3, 10}, Frame -> True, Axes -> False, 
      PlotLabel -> "Potential Barrier", 
      FrameLabel -> {"x", "ϵ"}, RotateLabel -> False]

Case 1: Energy particle E < V

We make the variable change ϵ ->  ((2 m ϵ)/(ℏ^2)). 
Lets deal with a particle caming from - ∞ to ∞. The proposed solution is:

\(ψ1 = \ A1\ Exp[\ I\ k1\ x] + B1\ Exp[\(-I\)\ k1\ x]; 
	

    ψ2 = \ A2\ Exp[\ \ α\ x] + B2\ Exp[\ \(-\ α\)\ x]; 

    		
ψ3 = \ \ A3\ Exp[\ I\ k1\ x]; 


    kToϵ\  = 
      \ {\ k1\  -> \ \@ϵ, α -> \ \@\(V2 - ϵ\)}; \)

Warning: It is a common mistake to set B2 = 0 at first sight

Parameters  =  {L -> 1, V2  -> 4 , ϵ -> 3 , A1 -> 1};

Now we apply the continuity and devivability conditions

\(\(Sol = 
      Solve[
	{ψ1\  == ψ2\  /. x -> 0\ , \ \ 
          D[ψ1, x] == D[ψ2, x]\  /. x -> 0\ , \ 
	
          \(ψ2\  == ψ3\  /. x -> L\)\  /. Parameters, \ \ 
          \(D[ψ2, x] == D[ψ3, x] /. x -> L\) /. Parameters}
, {A2, 
          B1, B2, A3}]; \)\)

\(ψ1Sol = \ 
      \(\(A1\ Exp[\ I\ k1\ x] + B1\ Exp[\(-I\)\ k1\ x] /. Sol\) /. 
          kToϵ\) /. Parameters; 
	

    ψ2Sol = \ 
      \(\(A2\ Exp[\ \ α\ x] + B2\ Exp[\ \(-\ α\)\ x] /. Sol
            \) /. kToϵ\) /. Parameters; 
		

    ψ3Sol = \ \ 
      \(\(A3\ Exp[\ I\ k1\ x] /. Sol\) /. kToϵ\) /. Parameters; \)

\(Show[
	
      Plot[\ Re[ψ1Sol]\ , {x, \(-10\), 0}, DisplayFunction -> Identity\ ]
        
, 
		
      Plot[\ Re[ψ2Sol]\ , {x, 0, L /. Parameters}\ , 
        DisplayFunction -> Identity]
, 
		
      Plot[\ Re[ψ3Sol]\ , {x, L /. Parameters, 10}\ , 
        DisplayFunction -> Identity]
, 
	
      Plot[\ ϵ\  /. Parameters, {t, \(-10\), 10}, Frame -> True, 
        PlotStyle -> {Dashing[{0.01, 0.01}]}, DisplayFunction -> Identity]
, 
      
	Plot[\ Evaluate[\(V[L, V2]\)[x] /. Parameters]\ , {x, \(-10\), 10}, 
        Frame -> True, Axes -> False, PlotLabel -> "Potential Barrier", 
        FrameLabel -> {"x", "ϵ"}, RotateLabel -> False, 
        DisplayFunction -> Identity]
, Frame -> True, 
      DisplayFunction -> $DisplayFunction, PlotLabel -> "Re[ψ]", 
      FrameLabel -> {"x", "ϵ"}, RotateLabel -> False, 
      Axes -> False]\)

Differential Operators

Hermite Polynomials

Hermite Polynomials can be obtained by nested application of a differential operator.

nM = 5; 

opH = 2*x*#1 - D[#1, x] & 

Expand[  NestList[opH, 1, nM]  ]

In order to be sure of the result let us use the Hermite Polynomials pre defined as HermiteH

Table[HermiteH[n,  x], {n,  0,  nM}]

Commutators

\(PO\  = \ 1/I\ \ D[#, x]&; 
XO\  = \ \ x\ #\ &; 

    \(Commutator[x]\)[XO, PO]\)

\(Lx\  = \ y\ D[#, z]\  - \ z\ D[#, y]\ &; 

    Ly\  = \ z\ D[#, x]\  - \ x\ D[#, z]&\ ; 


    \(Commutator[{x, y}]\)[Lx, Ly]\)

Angular Momentum.

Definition of Angular Momentum Operators

Lets take the special case of a system of one particle of spin 1/2 . As usual we take a Basis of both the L^2 and the component z of L.
When SetAngularMomentumOperators is called other operators like JPLUS and JMINUS are defined.

\(SetAngularMomentumOperators[L]
\)SetAngularMomentumKets[L, 1/2]

J2[L][ Ket[L, {1/2, 1/2}] ]\(\(Jz[L]\)[\ Ket[L, {1/2, 1/2}]\ ]
\)JPLUS[L][ Ket[L, {1/2, -1/2}] ]\(\(JMINUS[L]\)[\ Ket[L, {1/2, \(-1\)/2}]\ ]
\)

Matrix Representation

\(HIP[
	Bra[L, {1/2, 1/2}], \ Ket[L, {1/2, \(-1\)/2}]
]\)

\(Table[
	

      HIP[
			Bra[L, {1/2, r}]\ , \ \(JPLUS[L]\)[\ Ket[L, {1/2, s}]]\ 

        		], 
	
	{r, \(-1\)/2, 1/2}, {s, \(-1\)/2, 1/2}]\)

\(\(Jx[L]\)[w_]\  = \ 
      1/2 \((\(JPLUS[L]\)[w]\  + \ \(JMINUS[L]\)[w])\)\  // Expand; 

    \(Jy[L]\)[w_]\  = \ 
      1/\((2  I)\) \((\(JPLUS[L]\)[w]\  - \ \(JMINUS[L]\)[w])\) // Expand; 
    \)

\(σx\  = 
      2\ Table[
		
          HIP[	
Bra[L, {1/2, r}], \ \(Jx[L]\)[\ Ket[L, {1/2, s}]]
			]

          	, {r, \(-1\)/2, 1/2}, {s, \(-1\)/2, 1/2}]\)

\(σy\  = 
      2\ Table[
			
          HIP[
Bra[L, {1/2, r}], \ \(Jy[L]\)[\ Ket[L, {1/2, s}]]
			]

          			, {r, \(-1\)/2, 1/2}, {s, \(-1\)/2, 1/2}]\)

\(σz\  = 
      2\ Table[
		
          HIP[	
Bra[L, {1/2, r}], \ \(Jz[L]\)[\ Ket[L, {1/2, s}]]

            			\ \ \ ]
			, {r, \(-1\)/2, 1/2}, {s, \(-1\)/2, 1/2}]\)

Pauli Matrix Properies

Commutator[σx, σy]Commutator[σy, σz]\(Commutator[σz, σx]
\)σx.σx\(σy.σy
\)σx.σy  + σy.σxσx.σz  + σz.σx

We know the  Lz  Basis, but we also would like to know the Lx or Ly Basis. Lets find the Ly Basis.

{Ket[{L, Ly}, -1/2], Ket[{L, Ly}, 1/2]} = 
      Eigenvectors[1/2 σy].{Ket[L, {1/2, -1/2}], 
          Ket[L, {1/2, 1/2}]};

Ket[{L, Ly}, -1/2] // KetNormalizationKet[{L, Ly}, +1/2] // KetNormalization

Angular Momentum Adding.

Orbital +  Electron Spin.

SetAngularMomentumOperators[L]\(SetAngularMomentumOperators[S]
\)\(\(SetAngularMomentumKets[L, #]&\)\ /@\ {0, 1, 2, 3}; 

    SetAngularMomentumKets[S, 1/2]; \)

JPLUS[S][ Ket[S, {1/2, -1/2}] Ket[L, {0, 0}] ]

\(\(\(JMINUS[L + S] = \ \((\(JMINUS[S]\)[#] + \ \(JMINUS[L]\)[#])\)\ &; 
    \)

\)\)

\(\({\ Ket[J, {3/2, 3/2}]\ , Ket[J, {3/2, 1/2}], Ket[J, {3/2, \(-1\)/2}], 
        Ket[J, {3/2, \(-3\)/2}]} = 

      NestList[
	JMINUS[L + S], \ Ket[S, {1/2, 1/2}]\ Ket[L, {1, 1}], 3]; 
    \)\)

\(\(\ Ket[J, {3/2, 1/2}]\  // KetNormalization
\)\)

Two Electrons.

SetAngularMomentumOperators[S1]\(SetAngularMomentumOperators[S2]
\)\(SetAngularMomentumKets[S1, 1/2]; 
SetAngularMomentumKets[S2, 1/2]; 
    \)

For      STot = 1

\(\(Ket[STot, {1, 1}]\  = \ Ket[S1, {1/2, 1/2}] Ket[S2, {1/2, 1/2}]; 


    JMINUS[S1 + S2]\  = \ \(JMINUS[S1]\)[#] + \(JMINUS[S2]\)[#]&; 


    
{Ket[STot, {1, 1}], Ket[STot, {1, 0}], Ket[STot, {1, \(-1\)}]} = 

      NestList[	JMINUS[S1 + S2], \ Ket[STot, {1, 1}]\ , 2]\ ; 


    Ket[STot, {1, 0}]\  // KetNormalization
\)\)

ClebschGordan Coefficients

There is another aproach to get the properly combination and it is by using the  ClebschGordan[{(Subscript["\!\(\*StyleBox[\"\\\"j\\\"\", \"TI\"]\)", 1]), (Subscript["\!\(\*StyleBox[\"\\\"m\\\"\", \"TI\"]\)", 1])}, {(Subscript["\!\(\*StyleBox[\"\\\"j\\\"\", \"TI\"]\)", 2]), (Subscript["\!\(\*StyleBox[\"\\\"m\\\"\", \"TI\"]\)", 2])}, {j, m}]  coefficients. Lets take an example.

\(Sum[
	\ 
      Ket[S1, {1/2, m1}] Ket[S2, {1/2, m2}]*
        ClebschGordan[{1/2, m1}, {1/2, m2}, {1, 0}]
	

      	, {m1, \(-1\)/2, 1/2}, {m2, \(-1\)/2, 1/2}]\)

Note.- Dont worry about the warnings

Angular Differential Operators

Diferential operators can be defined naturally as pure functions.

RO  =  (x e[x] + y e[y] + z e[z]) #  &\(PO\  = \(1/I\ Grad[#]&\)\ \  /. \ Grad[Cartesian]
\)\(P2O = \ \ \(PE[{x, y, z}]\)[PO, PO]
\)LO  =  PV[{x, y, z}][RO, PO]

We can have these operators in other curvilinear systems

\(PS\  = \ \(1/I\ Grad[#]&\)\ \  /. \ Grad[Spherical]
	\)\(PS2\  = \ \(PE[{r, th, ph}]\)[PS, PS]
\)

RS  =  r e[r] # &LS =  PV[{r, th, ph}][RS, PS]

PE[{r, th, ph}][LS, LS]

Podemos extraer la componente x de la siguiente manera

\(\(exSp\  = \ \(e[x] #\ &\) /. \ TransformVersors[Cartesian, Spherical]; 
    

LX\  = \(PE[{r, th, ph}]\)[\ \ exSp\ , LS\ ]\ \)\)

The conmutator is also valid for vectorial differential operators  (Not Fully Tested)

Commutator[{x, y, z}][RO, PO]

Angular Momentum in Spherical Coordinates

The spherical coordinates are very special and also very important representation of  the Angular Momentum. Lets define the Angular Momentum Ortonormal Basis

\(SetAngularMomentumOperators[L1]
\)SetAngularMomentumKets[L1, 1]

We would need some graphical functions

<< Graphics`ParametricPlot3D`

Lets find the d orbital with  m=0 and get its plot

G[d0]  =  
      Abs[ HIP[Bra[Spherical, {th, ph}], Ket[L1, {2, 0}]] ]  // Evaluate

\(PointParametricPlot3D[\ 
	
      G[d0] {Sin[th] Cos[ph], Sin[th] Sin[ph], Cos[th]}	, {th, 0, Pi}, {ph, 
        0, 2  Pi}, PlotPoints -> 40\ \ , Boxed -> False]\)

Simple Harmonic Oscillator

Annihilation Creation Operators

First we define the discrete ortonormal Basis

SetDiscreteOB[ Ket[HO, {0, 1}]]

Now is time to define the Annihilation and Creation operators

\(Annihilation[\ Ket[HO, m_]] := \ Sqrt[m] Ket[HO, m - 1]; 

    Annihilation[Ket[HO, 0]] = 0; 


    Creation[\ Ket[HO, m_]] := \ Sqrt[m + 1] Ket[HO, m + 1]; 

    Creation[0] = 0; 

SetLinearOperator[Creation, Ket]\)SetLinearOperator[Annihilation, Ket]

Annihilation[a (Ket[HO, k]/b + Ket[HO, m]/c)]

\(Table[
	HIP[
	Bra[HO, i], Creation[\ Ket[HO, j]\ ]
		\ \ \ \ ]

      	\ \ \ \ \ \ , {i, 0, 4}, {j, 0, 4}]
\)\(Table[
	
      HIP[
	Bra[HO, i], Annihilation[\ Ket[HO, j]\ ]
		\ \ \ ]

      	\ \ \ \ \ \ , {i, 0, 4}, {j, 0, 4}]\)

\(NPhonones\  = \ Creation[Annihilation[#]]&; 


    Table[
\ \ 
      HIP[	
\ \ \ Bra[HO, i], NPhonones[Ket[HO, j]]	\ 
	\ \ \ \ \ \ \ ]

      	\ \ \ \ \ \ , {i, 0, 4}, {j, 0, 4}]\)

Defining the position and the momentum operators

\(X\  = \ 1/Sqrt[2] \((Annihilation[#]\  + \ Creation[#]\ \ )\)&; 

    P\  = \ 1/\((I\ Sqrt[2])\) \((Annihilation[#]\  - \ Creation[#])\)&; 		
    

HIP[
	Bra[HO, k], X[X[\ Ket[HO, k]\ ]]
]\  // Expand
\)\(HIP[
	Bra[HO, k], P[P[\ Ket[HO, k]\ ]]
] // Expand\ \ \)

Now a little more tedious  calculus by hand, but not by computer!

\(X6\  = \ Composition[X, X, X, X, X, X]; 
X6[Ket[HO, k]] // Expand\)

\(HIP[
	Bra[HO, k], \ \ X6[\ Ket[HO, k]]
]\  // Expand\)

Defining the Hamiltonian of the Simple Harmonic Oscilator

\(HamHO\  = \ NPhonones[#]\  + \ 1/2  #\ &; 


    HamHO[Ket[HO\ , n]]\  // Factor\)

Wave function

In order to see their wave function we can use the predefined ket :
Ket[ HarmonicOscillator  ,{n}]

\(HIP[
	Bra[\ Cartesian, {x}], Ket[\ HarmonicOscillator\ \ , {n}]\ 

      	\ \ \ ]\)

\(HOWF\  = 
      Table[
	\ 
        HIP[Bra[\ Cartesian, {x}], Ket[\ HarmonicOscillator\ \ , {n}]\ ]

        	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ , {n, 0, 5}]\)

Plot[Evaluate[HOWF], {x, -5, 5}, Frame -> True, 
      FrameLabel -> {"x", "ψ"}, 
      PlotLabel -> "Wave functions of the Harmonic Oscilator "]

Hidrogen Atom Without Spin

We can obtain the Hidrogen Atom Wave Function in Spherical Coordinates by the use of predefined Ket  " Ket[HidrogenWS,{n,l,m}]  "

\(HIP[
Bra[Spherical, {r, th, ph}], Ket[HidrogenWS, {3, 2, 0}]] // 
      HoldForm\)

\(Hidrod0\  = \ 
      \(r^2\ HIP[
Bra[Spherical, {r, th, ph}], Ket[HidrogenWS, {3, 2, 0}]

              \ \ \ \ \ \ \ ] /. \ 
          TransformVariables[Spherical, Cartesian]\) /. x -> 0\)

Now we can show this probability function cut along the y-z plane

DensityPlot[Abs[Hidrod0]^2, {y, -15, 15}, {z, -15, 15}, 
      Mesh -> False, PlotPoints -> 50]

Time Independient Perturbation Theory

Nondegenerate Case

Implementation

We are going to apply the particulary method that can be see in "J.J. Sakurai, Modern Quantum Mechanics, Addison Wesley Publishing Company Inc., 1985  " , because it is more suitable for symbolic computation as we'll see below.

The Hamiltonian to be solved is

((Subscript[H, 0]) + λ V) | n >=  (Subscript[E, n]) | n >

The system of equations to be solved given |(n^((0))) >   are

\(\(\(|\(n > \)\  = \ |\) n\^\((0)\) > 
        \(+\ \ \(ϕ\_n\/\(E\_n\^\((0)\) - H\_0\)\)\) 
          \((λ\ V\  - Δ\_n)\) | \(n > \)
		\)\)

\(\(	Δ\_n = λ < n\^\((0)\) | \(n > \)\)\)

Which (Subscript[ϕ, n]) stands for the Complementary Projection Operator. Now lets put these equation in computational code.

\(\(SetDiscreteOB[\ 	
      Ket[a[0], {0, 1}]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 


    SetLinearOperator[CPO, Ket]; 

SetLinearOperator[H0, Ket]; 


    SetLinearOperator[V, Ket]; 
\)\)

\(\(CPO[\ \ Ket[a[0], _\ ]\ \ ]\  = 0; 


    Format[CPO[w_]]\  = \ \ 
      SequenceForm[\ 
        ϕ/\((\ Subscript["E", n] - \ Subscript[H, 0])\)\ , w]; 

    \)\)

Then the equation (2 ) becomes to be

\(KetEquation\  = \ 
      Ket[a, n] == 	
        Ket[\ a[0]\ , n]\  + \ \ 

          CPO[\(\((λ\ V[#]\  - \ Δ[n] #)\)&\)[
              Ket[a, n]]]\)

And the equation (3)

ΔEnergy  =  
      Δ[n] ==  
        λ HIP[   Bra[a[0], n]  , V[ Ket[a, n] ] ]

Now we need to expand both  | n >   and   Δ[n]  and to do so we define rules of expansion.

ExpandΔ[n][m_]  =  
      Δ[n] ->  
        Sum[ (λ^i) Δ[i] , {i, 1, m}]

ExpandDelta[m_]  =  
      Delta[n] ->  Sum[ Lambda^i  Delta[i] , {i, 1, m}]

ExpandKet[a, n][m_]  =  
      Ket[a, n] ->  Sum[ (λ^i) Ket[a[i], n] , {i, 0, m}]

ExpandKet[m_]  =  
      Ket[a, n] ->  Sum[ Lambda^i Ket[a[i], n] , {i, 0, m}]

Note.- "m"   must be at least the value of the  maximun order required.

Applying these rules over Δ[n]

\(\(ΔEnergyExpansion\  = \ 
      \(ΔEnergy\  /. \ \(ExpandΔ[n]\)[3]\) /. 
        \(ExpandKet[a, n]\)[3]\ \ 
\)\)

Right below the function that give us the (λ-th) coefficient

\(TerminoOrden[w_, {λ_^n_.}] := 
      Module[{}, 

          Coefficient[
			\ Apply[Subtract, w // ExpandAll], λ^n]]
           == 0\)

And we certainly have the values for the n-th shitf energy

TerminoOrden[ΔEnergyExpansion, {λ}]TerminoOrden[ΔEnergyExpansion, {λ^2}]TerminoOrden[ΔEnergyExpansion, {λ^3}]

Applying the expansion rules over  the ket equation (2)

KetEquationExpansion  =  
      KetEquation /.  ExpandΔ[n][4] /. 
        ExpandKet[a, n][4]  ;

TerminoOrden[KetEquationExpansion, {λ}]  TerminoOrden[KetEquationExpansion, {λ^2}]  \(TerminoOrden[KetEquationExpansion, {λ^3}]\ \ 

\)

Now we would like a more compact form to display the  n-th order ket and shift energy correction.

PK =PerturbationKets[Generic[a, V, CPO, Delta, n], 4]  // TableForm

PerturbationShiftEnergy[Generic[a,  V,  Delta,  n],  3]  // 
        TableForm