I have several questions concerning manipulations with creation and annihilation operators. I would be grateful if besides answering my questions you would give me some suggestions for improving the code.
Background of the physic problem
I have a crystal lattice with several sites per unit cell. Every site can holds one, two or none of electrons. Electron at every site can be at one of two states: spin up and spin down. Electrons are Fermi particles so anticommutative basis is used. Index 'i' (1<=i<=f) corresponds to i's lattice cell, index 'j' (1<=j<=3) corresponds to j's site within i'th lattice cell, index 'sigma' corresponds to spin of the electron (if sigma=1 then spin is down, if sigma=2 then spin is up).
>with Physics :
>Setup mathematicalnotation = true :
>Physics:-Setup anticommutativeprefix = Psi :
>#Below the annihilation operators are defined. They depend on indicies [i,j,sigma].
>am:=proc(i, j, sigma)
>local k: : posint :
>Annihilation(Psi, k, notation = explicit);
>#Below the creation operators are defined. They depend onindicies [i, j, sigma].
>ap:=proc(i, j, sigma)
>local k :: posint :
>Creation(Psi, k, notation = explicit);
>#Here the occupation number at state i, j, sigma is programmed
>nn:=proc( i, j, sigma)
>ap(i, j, sigma).am(i, j, sigma)
>#Now I want to define operator of total number of particles.
># I have tried the:
>#N:=sum( sum( sum( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;
>#N:=Sum( Sum( Sum( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;
>#but both give an error.
>#The only working way is
>f:=3; # 'f' have to be defined explicitly
>N:=add( add( add( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;
Is there a way to define operator 'N' with indefinite 'f'?
Is there a way to have operators nn in output of N rather then a^+a?
Is there a tool for simplification of the operators?
I'd like to have something like