Implementation of power operation between a qnumber and a rnumber

[labay11]

In driven-dissipative bosonic systems it is possible to have high nonlinear terms in the Hamiltonian and in the dissipative parts. For instance, n-photon driving and disspation. However, the package does not allow to perform operations of the type a^n where a is a qnumber representing the annihilation operator in Fock space and n a rnumber.

[david-pl]

While I understand the motivation, this would be very difficult to implement and even then it would only work for some specific cases. I general, we could allow expressions of the form a^n, but they would only work up to a certain point. The problem is that the cumulant expansion is dependent on the order of an average. For example, the second order expansion of the average $\langle a^n\rangle$ changes depending on the value of $n$. This means that when you derive a set of equations for such a system, the value of $n$ changes the number of equations required to complete the set of equations and solve the system. So the latest point at which you have to choose a concrete value for n is before you apply the cumulant expansion and complete the set of equations. Quite frankly, you might as well just choose a specific value right away and re-run the program for a different value when required.