Abstract
Using the positive P-representation of the density matrix, we develop an algorithm for calculating the quantum many-body nonlinear response functions of a system of bosons driven impulsively by external fields. The formalism maps the quantum time evolution of N boson degrees of freedom into a stochastic dynamics of 4N classical degrees of freedom. The first- and the third-order response functions are calculated by propagating the parameters of the P-representation using a set of coupled Langevin equations with multiplicative noise. These parameters serve as classical variables. Two classical ways for computing the response functions are presented. In the nonequilibrium method, an observable is calculated for weak impulsive pulses, and the response functions are obtained by taking its derivatives with respect to the pulse amplitudes. In the alternative, equilibrium simulation, the response functions are expressed in terms of time-correlation functions involving the P-representation parameters and stability matrices representing the perturbation of the trajectories. The stability matrices can be propagated simultaneously with the Langevin equations for the parameters. The formalism is generalized for a many-body boson system coupled to a harmonic bath.