# Second Quantization in Matlab

This question may be more suited for physics.stackexchange, but I saw this post was recommended for StackOverflow or Computational Science, so I'm asking my question here.

I am trying to write a program in Matlab to evaluate an expression with fermionic annihilation and creation operators. Define $$\{f, g \} \equiv fg + gf$$ then the fermionic second quantization operators satisfy $$\{a_i,a_j \} = \{ a_i^\dagger, a_j^\dagger\} = 0 \qquad \{a_i^\dagger, a_j\} = \delta_{ij}$$ $\forall i, j \in \{1, 2, ..., N \}$ for some given $N$.

$W(t)$ is a $2N \times 2N$ matrix that I know that gives me the time evolution of the operators in the Heisenberg picture; \begin{align*} a_i(t) &= \sum_{j=1}^N \left(W_{2i-1,2j-1}(t)a_j + W_{2i-1,2j}(t)a_j^\dagger\right)\\ a_i^\dagger(t) &= \sum_{j=1}^N \left(W_{2i,2j-1}(t)a_j + W_{2i,2j}(t)a_j^\dagger\right) \end{align*} i.e. $$\left(\begin{matrix}a_1(t)\\a_1^\dagger(t)\\\vdots\\a_N(t)\\a_N^\dagger(t) \end{matrix}\right) = W(t)\left(\begin{matrix}a_1\\a_1^\dagger\\\vdots\\a_N\\a_N^\dagger \end{matrix}\right)$$

I want to write a function that takes in $W(t)$ (evaluated at a given $t$) and $j$ and returns $$\langle \psi | \left(\prod_{i=1}^{j-1}\left(1-2a_i^\dagger(t)a_i(t) \right) \right)\left(a_j(t)^\dagger + a_j(t) \right) | \psi \rangle$$ where $$| \psi \rangle = \left(\prod_{i=1}^N\frac{1+a_i^\dagger}{\sqrt{2}} \right)|0\rangle = 2^{-N/2}\left(|0...0\rangle + |0...1 \rangle + |0...1 \ 0 \rangle + |0 ...1 \ 1 \rangle + ... + |1...1 \rangle \right)$$ $|\psi\rangle$ is just a normalized, equally weighted linear combination of all states. Note that it does not change with time. The time evolution is encoded in the operators.

My thought was to use Matlab's symbolic algebra, and try and code for what each operator does to a state. But the symbolic algebra does not always preserve order, which is crucial because the operators don't always commute. Also, trying to just naively go through is so computationally expensive and grows exponentially with $N$.

Any ideas?

• For the non-commutative part and computer algebra, you can use modules/packages that handle non-commutative algebras. On the other hand, if you think that numerically it can become intractable, because it grows exponentially, the same would happen symbolically. And symbolic calculations are more expensive, normally. – nicoguaro Jul 18 '17 at 13:43
• Can you suggest a module/package to use? – Joe Jul 18 '17 at 13:59
• Not really, I used Maple for Dirac's matrices (that form a Clifford algebra) a long time ago. I know that Maxima has the atensor package that allows manipulation of noncommutative algebras. Maybe Cadabra has something ... It might be a better idea to check the Comparison of CAS in Wikipedia. – nicoguaro Jul 18 '17 at 19:34