  We present a method of calculating the interacting S-matrix to an arbitrary
perturbative order for a large class of boson interaction Lagrangians. The
method takes advantage of a previously unexplored link between the $n$-point
Green's function and a certain system of linear Diophantine equations. By
finding all nonnegative solutions of the system, the task of perturbatively
expanding an interacting $S$-matrix becomes elementary for any number of
interacting fields, to an arbitrary perturbative order (irrespective of whether
it makes physical sense) and for a large class of scalar boson theories. The
method does not rely on the position-based Feynman diagrams and promises to be
extended to many perturbative models typically studied in quantum field theory.
Aside from interaction field calculations we showcase our approach by expanding
a pair of Unruh-DeWitt detectors coupled to Minkowski vacuum to an arbitrary
perturbative order in the coupling constant. We also link our result to Hafnian
as introduced by Caianiello and present a method to list all (2n-1)!! perfect
matchings of a complete graph on 2n vertices.
