  Correlation functions in the O(n) models below the critical temperature are
considered. Based on Monte Carlo (MC) data, we confirm the fact stated earlier
by Engels and Vogt, that the transverse two-plane correlation function of the
O(4) model for lattice sizes about L=120 and small external fields h is very
well described by a Gaussian approximation. However, we show that fits of not
lower quality are provided by certain non-Gaussian approximation. We have also
tested larger lattice sizes, up to L=512. The Fourier-transformed transverse
and longitudinal two-point correlation functions have Goldstone mode
singularities in the thermodynamic limit at k --> 0 and h=+0, i.e., G_perp(k) =
a k^{-lambda_perp} and G_parallel(k) = b k^{-lambda_parallel}, respectively.
Here a and b are the amplitudes, k is the magnitude of the wave vector. The
exponents lambda_perp, lambda_parallel and the ratio b M^2/a^2, where M is the
spontaneous magnetization, are universal according to the GFD (grouping of
Feynman diagrams) approach. Here we find that the universality follows also
from the standard (Gaussian) theory, yielding b M^2/a^2 = (n-1)/16. Our MC
estimates of this ratio are 0.06 +/- 0.01 for n=2, 0.17 +/- 0.01 for n=4 and
0.498 +/- 0.010 for n=10. According to these and our earlier MC results, the
asymptotic behavior and Goldstone mode singularities are not exactly described
by the standard theory. This is expected from the GFD theory. We have found
appropriate analytic approximations for G_perp(k) and G_parallel(k), well
fitting the simulation data for small k. We have used them to test the
Patashinski--Pokrovski relation and have found that it holds approximately.

Statistical Mechanics