  We rigorously derive weakly nonlinear relation between cosmic density and
velocity fields up to third order in perturbation theory. The density field is
described by the mass density contrast, $\de$. The velocity field is described
by the variable $\te$ proportional to the velocity divergence, $\te = -
f(\Omega)^{-1} H_0^{-1} \nabla\cdot\bfv$, where $f(\Omega) \simeq
\Omega^{0.6}$, $\Omega$ is the cosmological density parameter and $H_0$ is the
Hubble constant. Our calculations show that mean $\de$ given $\te$ is a third
order polynomial in $\te$, $\lan \de \ran|_{\te} = a_1 \te + a_2 (\te^2 -
\s_\te^2) + a_3 \te^3$. This result constitutes an extension of the formula
$\lan \de \ran|_{\te} = \te + a_2 (\te^2 - \s_\te^2)$, found by
Bernardeau~(1992) which involved second order perturbative solutions. Third
order perturbative corrections introduce the cubic term. They also, however,
cause the coefficient $a_1$ to depart from unity, in contrast with the linear
theory prediction. We compute the values of the coefficients $a_p$ for
scale-free power spectra, as well as for standard CDM, for Gaussian smoothing.
The coefficients obey a hierarchy $a_3 \ll a_2 \ll a_1$, meaning that the
perturbative series converges very fast. Their dependence on $\Omega$ is
expected to be very weak. The values of the coefficients for CDM spectrum are
in qualitative agreement with the results of N-body simulations by Ganon et al.
(1996). The results provide a method for breaking the $\Omega$-bias degeneracy
in comparisons of cosmic density and velocity fields such as IRAS-POTENT.
