Correlation functions in the $\mathcal{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 $\mathcal{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_⊥$(k)$≃ak^{−λ_⊥}$ and $G_{||}($k$)≃bk^{−λ_{||}}$, respectively. Here $a$ and $b$ are the amplitudes, $k$=|k| is the magnitude of the wave vector k. The exponents $λ_⊥$, $λ_{||}$ and the ratio $bM^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 $bM^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_⊥$(k) and $G_{||}$(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.