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 <PERSON> 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. 281
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? 82