I'm considering these functions $f,g,h,j$, which all share the same form, of the boson creation and annihilation operators:
$$ f(\hat{a}^{\dagger} ,\hat{a}) = \sum_{n,m} c_{n,m}(\hat{a}^{\dagger})^n (\hat{a})^m \hspace{3cm} [\hat{a} , \hat{a}^{\dagger}] = 1$$
and trying to calculate the following commutator:
$$[f(\hat{a}^{\dagger} ,\hat{a}) \hat{a} g(\hat{a}^{\dagger} ,\hat{a}), h(\hat{a}^{\dagger} ,\hat{a}) \hat{a}^{\dagger} j(\hat{a}^{\dagger} ,\hat{a})]$$
$$= \sum_{n,m, \alpha, \beta, i, j, \rho, \sigma} c_{n,m} c_{ \alpha, \beta} c_{ i, j} c_{ \rho, \sigma} [ (\hat{a}^{\dagger})^n (\hat{a})^m \hat{a} (\hat{a}^{\dagger})^{\alpha} (\hat{a})^{\beta} , (\hat{a}^{\dagger})^i (\hat{a})^j \hat{a}^{\dagger} (\hat{a}^{\dagger})^{\rho} (\hat{a})^{\sigma} ]$$
$$= \sum_{I} C_I [ (\hat{a}^{\dagger})^n (\hat{a})^{m+1} (\hat{a}^{\dagger})^{\alpha} (\hat{a})^{\beta} , (\hat{a}^{\dagger})^i (\hat{a})^j(\hat{a}^{\dagger})^{\rho+1} (\hat{a})^{\sigma} ]$$
where $I := \{ n,m, \alpha, \beta, i, j, \rho, \sigma \}$ is an index set and $C_I := c_{n,m} c_{ \alpha, \beta} c_{ i, j} c_{ \rho, \sigma} $ is the amalgamation of all the constants.
In another post I've seen the result:
$$[\hat{a}^n,(\hat{a}^{\dagger})^m]~=~\sum_{k=1}^{\min(n,m)} \frac{n!m!}{(n-k)!(m-k)! k!}(\hat{a}^{\dagger})^{m-k}\hat{a}^{n-k},$$
which simplifies the calculations substantially, however, it still needs to be used a tally of $24$ times.
It's very big and I wonder if there is a general result about these commutators that allows me to compute that commutator even faster.
