Hiking a generalized Dyck path: A tractable way of calculating multimode boson evolution operators

A time evolution operator in the interaction picture is given by exponentiating an interaction Hamiltonian $H$. Important examples of Hamiltonians, often encountered in quantum optics, condensed matter and high energy physics, are of a general form $H=r(A^\dagger-A)$, where $A$ is a multimode boson operator and $r$ is the coupling constant. If no simple factorization formula for the evolution operator exists, the calculation of the evolution operator is a notoriously difficult problem. In this case the only available option may be to Taylor expand the operator in $r$ and act on a state of interest $\psi$. But this brute-force method quickly hits the complexity barrier since the number of evaluated expressions increases exponentially. We relate a combinatorial structure called Dyck paths to the action of a boson word (monomial) and a large class of monomial sums on a quantum state $\psi$. This allows us to cross the exponential gap and make the problem of a boson unitary operator evaluation computationally tractable by achieving polynomial-time complexity for an extensive family of physically interesting multimode Hamiltonians. We further test our method on a cubic boson Hamiltonian whose Taylor series is known to diverge for all nonzero values of the coupling constant and an analytic continuation via a Pad\'e approximant must be performed.