Implement operator-based PT

Right now defining the perturbative problem requires specifying the basis of the Hilbert space, which doesn't apply naturally to the second-quantized problems. This means that e.g. going from a Hubbard model to Heisenberg, computing ring exchange, or more broadly dealing with big interacting Hamiltonians requires more work than necessary, and needs the users deal with the concepts that they don't usually apply.

Below I outline a rough and incomplete idea of what we can do better.

As an alternative, we can formulate the entire perturbative problem in the language of operators. Roughly this means that we decompose every single operator into a weighed sum of normally ordered terms of the form O_i = \left[\prod_{i=1}^{N_F} (a^\dagger_i)^{n_{i+}}(a_i)^{n_{i-}}\right]\times\left[ \prod_{i=1}^{N_B} (b^\dagger_i)^{m_{i+}}(b_i)^{m_{i-}}\right] \times \left[\prod_{i=1}^{n_\sigma} \sigma_i^{(\alpha_i)}\right], where n_{i\pm} \in \{0, 1\}, m_{i\pm} \in \mathbb{N}, and \alpha_i \in \{0, x, y, z\}.

In terms of code we can store this information in a dictionary where keys are all the n_i, m_i, and \alpha_i and the values are the amplitudes of the Hamiltonian terms.

Implementing addition and multiplication of such operators is relatively straightforward, which only leaves open a question of how we would solve the analog of the Sylvester's equation, which now amounts to finding the least squares solution to [H_0, V_n] = Y_n. This, I think is only possible if H_0 is a sum of commuting operators acting on different degrees of freedom. If that condition is fulfilled, then the Sylvester's equation does not mix operators containing different sets of second-quantized operators. This means that if Y_n has a finite number of terms x_i O_i, then V_n will also always have a finite number of terms, each containing products of operators that are acting on the same fermions, bosons, and spins, as those appearing in Y_n.