  We show that the formation of large-scale structures through gravitational
instability in the expanding universe can be fully described through a
path-integral formalism. We derive the action S[f] which gives the statistical
weight associated with any phase-space distribution function f(x,p,t). This
action S describes both the average over the Gaussian initial conditions and
the Vlasov-Poisson dynamics. Next, applying a standard method borrowed from
field theory we generalize our problem to an N-field system and we look for an
expansion over powers of 1/N. We describe three such methods and we derive the
corresponding equations of motion at the lowest non-trivial order for the case
of gravitational clustering. This yields a set of non-linear equations for the
mean $\fb$ and the two-point correlation G of the phase-space distribution f,
as well as for the response function R. These systematic schemes match the
usual perturbative expansion on quasi-linear scales but should also be able to
handle the non-linear regime. Our approach can also be extended to non-Gaussian
initial conditions and may serve as a basis for other tools borrowed from field
theory.
