We investigate an alternative approach to the study of the solvent response to the sudden change in the charge distribution of a solute molecule. The theory avoids the assumption that the response induced in the solvent is linear with respect to the solute perturbation. Our method focuses on the nonequilibrium characteristic function gU (α;t) for the solvent contribution to the vertical energy gap, which provides a link between the averages measured in the solvation dynamics experiment and the molecular description of the nonequilibrium solvation process. We take advantage of the Kawasaki form of the nonequilibrium distribution function, which is valid only in the case of jump perturbations, to express the characteristic function as a ratio of two partition functions defined in terms of complex-valued, time-dependent, many-body effective Hamiltonians. We then apply the theory of the generalized Langevin equation to cast the partition functions in terms of approximate two-body additive effective Hamiltonians, in a way that enables us to exploit well known nonlinear equilibrium integral equation methodologies to investigate the process of nonequilibrium solvation on a molecular scale. To test the performance of our simplest approximation for gU (α;t) we report calculations of the nonequilibrium solvation time correlation function and of the evolution of the solvation structure for a model system (first studied by Fonseca and Ladanyi by molecular dynamics computer simulations) that displays important nonlinear effects.