  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.
