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†−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 ψ. 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  ψ. 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é approximant must be performed.