We develop a formalism for particle production in a field theory coupled to a strong time-dependent external source. An example of such a theory is the Color Glass Condensate. We derive a formula, in terms of cut vacuum-vacuum Feynman graphs, for the probability of producing a given number of particles. This formula is valid to all orders in the coupling constant. The distribution of multiplicities is non-Poissonian, even in the classical approximation. We investigate an alternative method of calculating the mean multiplicity. At leading order, the average multiplicity can be expressed in terms of retarded solutions of classical equations of motion. We demonstrate that the average multiplicity at next-to-leading order can be formulated as an initial value problem by solving equations of motion for small fluctuation fields with retarded boundary conditions. The variance of the distribution can be calculated in a similar fashion. Our formalism therefore provides a framework to compute from first principles particle production in proton-nucleus and nucleus-nucleus collisions beyond leading order in the coupling constant and to all orders in the source density. We also provide a transparent interpretation (in conventional field theory language) of the well known Abramovsky-Gribov-Kancheli (AGK) cancellations. Explicit connections are made between the framework for multi-particle production developed here and the framework of Reggeon field theory.