  The mass density distribution of Newtonian self-gravitating systems is
studied analytically in field theoretical method. Modeling the system as a
fluid in hydrostatical equilibrium, we apply Schwinger's functional derivative
on the average of the field equation of mass density, and obtain the field
equation of 2-point correlation function $\xi(r)$ of the mass density
fluctuation, which includes the next order of nonlinearity beyond the Gaussian
approximation. The 3-point correlation occurs hierarchically in the equation,
and is cut off by the Groth-Peebles anzats, making it closed. We perform
renormalization, and write the equation with three nonlinear coefficients. The
equation tells that $\xi $ depends on the point mass $m$ and the Jeans
wavelength scale $\lambda_{0}$, which are different for galaxies and clusters.
Applying to large scale structure, it predicts that the profile of $\xi_{cc} $
of clusters is similar to $\xi_{gg}$ of galaxies but with a higher amplitude,
and that the correlation length increases with the mean separation between
clusters, i.e, a scaling behavior $r_0\simeq 0.4d$. The solution yields the
galaxy correlation $\xi_{gg}(r) \simeq (r_0/r)^{1.7}$ valid only in a range
$1<r<10 \,h^{-1}$Mpc. At larger scales the solution $\xi_{gg} $ deviates below
the power law and goes to zero around $\sim 50 \, h^{-1}$Mpc, just as the
observations show. We also derive the field equation of 3-point correlation
function in Gaussian approximation and its analytical solution, for which the
Groth-Peebles ansatz with $Q= 1$ holds.

Astrophysics of Galaxies