Skip to main content

Symbolic programming

The beauty of Wolfram Language is its flexibillity. One can shape the language making a dialect of it for a specific problem to solve

Download original notebook

Wavefunctions​

Let us use built-in syntax for Ket and Bra vectors. Apart from the decorations, there is no predefined rules for those expressions

{Ket[1], Bra[1]}
{(*BB[*)(Ket[1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*),(*BB[*)(Bra[1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*)}

Let us define an ortogonallity rule on Dot operation aka

⟨Ms∣Ms′⟩=δMs,Ms′\langle Ms | Ms^\prime \rangle = \delta_{Ms, Ms^\prime} 
Ket /: (*BB[*)(Bra[Ms1_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*).(*BB[*)(Ket[Ms2_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := KroneckerDelta[Ms1, Ms2]
Ket /: Conjugate[(*BB[*)(Ket[Ms2_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)] :=  (*BB[*)(Bra[Ms2])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMaeiRJfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABBAEzo="*)(*]BB*)

now test it in the basis of a spin S=2S=2

∣S=2,Ms⟩| S=2, Ms \rangle 
basis = Table[(*BB[*)(Ket[Ms])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*), {Ms, -2,2}];

Table[Conjugate[i].j, {i, basis}, {j, basis}] 
{{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}}

or more appealing view

% // MatrixForm 
((*GB[*){{1(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)1(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)1(*|*),(*|*)0(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)1(*|*),(*|*)0}(*||*),(*||*){0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)0(*|*),(*|*)1}}(*]GB*))

Ladder operators​

Let us define a typical lowering and raising operators

Ket /: Sp.(*BB[*)(Ket[Ms_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := (*SqB[*)Sqrt[2(2+1) - (Ms+1)Ms](*]SqB*) (*BB[*)(Ket[Ms+1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
Ket /: Sm.(*BB[*)(Ket[Ms_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := (*SqB[*)Sqrt[2(2+1) - (Ms-1)Ms](*]SqB*) (*BB[*)(Ket[Ms-1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
Ket /: Sz.(*BB[*)(Ket[Ms_])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*) := Ms (*BB[*)(Ket[Ms])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)

Sx = (*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)(Sp + Sm);
Sy =  (*FB[*)((1)(*,*)/(*,*)(2I ))(*]FB*)(Sp - Sm);

Let's test

Sp.(*BB[*)(Ket[1])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*)
2 ((*BB[*)(Ket[2])(*,*)(*"1:eJxTTMoPSmNiYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFxTxAMe/UEpfU5PyixJL8IohKFiARVJqTGgxieKQmpgQzQ9QBABLAE1g="*)(*]BB*))

Linear algebra​

In Woflram Language a linear properties of Dot operator are not predefined, i.e.

a . (2 b) === 2 a.b 
a . (b + c) === a.b + a.c 
False
False

But it is relatively easy to define them

Unprotect[Dot];
Unprotect[Conjugate];

Conjugate[a_ b_] := Conjugate[a] Conjugate[b]
Conjugate[a_ + b_] := Conjugate[a] + Conjugate[b]

Dot[a_, b_?NumericQ] := b a
Dot[b_?NumericQ, a_] := b a
Dot[a_ + b_, c_] := a.c + b.c 
Dot[c_, a_ + b_] := c.a + c.b  
Dot[a_, b_?NumericQ c_] := b a.c
Dot[b_?NumericQ c_, a_] := b c.a

Protect[Dot];
Protect[Conjugate];

Now it is better

a . (2 b) === 2 a.b 
a . (b + c) === a.b + a.c 
True
True

Spin Hamiltonian​

For example, one can try to solve a well-known spin-hamiltonian for a spin system with an uniaxial anisotropy in magnetic field perpendicular to it

H^=S^xh−dS^z2\hat{H} = \hat{S}_x h - d \hat{S}_z^2

If d is negative, the spin system has easy-plane anisotropy, i.e. it is more preferable for a spin to lie in xy plane.

In a code it looks like this

H[h_, d_] := Sx h - d Sz.Sz

Let us find an average projection of a spin to the field direction

SxMatrix = Table[Conjugate[i].Sx.j, {i, basis}, {j,basis}];

Compute for each field the ground state and find an average of Sx

estimate[d_] := Table[
  With[{
    matrix = Table[Conjugate[i].H[h, d].j, {i, basis}, {j,basis}]
  },
    With[{
      system = Eigensystem[matrix]
    },
      With[{
        gState = {SortBy[system // Transpose, First] // First // Last}
      },
        {h, Conjugate[gState].SxMatrix.Transpose[gState] // First // First} 
      ]
    ]
  ]
, {h, 0, 5, 0.1}];


ListLinePlot[{estimate[-1], estimate[1]}, PlotLegends->{"easy-plane", "easy-axis"}, AxesLabel->{"h", "Sx"}]
(*VB[*)(FrontEndRef["f2fc3c79-a9f4-402a-9923-10790735cd7e"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKpxmlJRsnm1vqJlqmmeiaGBgl6lpaGhnrGhqYWxqYG5smp5inAgCEJhVP"*)(*]VB*)

One can modify it and estimate the projection for different anisotropy term d magnitudes

Taking it from 0 up to ~h

Table[
  With[{
    matrix = Table[Conjugate[i].H[h, d].j, {i, basis}, {j,basis}]
  },
    With[{
      system = Eigensystem[matrix]
    },
      With[{
        gState = {SortBy[system // Transpose, First] // First // Last}
      },
        {h, Conjugate[gState].SxMatrix.Transpose[gState] // First // First} 
      ]
    ]
  ]
, {d, 0, 10, 1}, {h, 0, 10, 0.2}] // ListLinePlot 
(*VB[*)(FrontEndRef["ea302655-c622-4044-a3b3-3cc4d7022ce3"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKpyYaGxiZmZrqJpsZGemaGJiY6CYaJxnrGicnm6SYGxgZJacaAwB4FhUn"*)(*]VB*)