If you have a harmonic oscillator in x, the ground state wave function is a Gaussian;
$$ H = {p^2\over 2} + {\omega^2 x^2\over 2} $$
$$ \psi_0(x) = e^{ - {\omega x^2\over 2}} $$
If you have two independent oscillators x,y;
$$ H = {p_x^2\over 2} + {p_y^2\over 2} + {\omega_1^2 x^2\over 2} + {\omega_2^2 y^2\over 2} $$
the ground state is a product:
$$ \psi_0(x,y) = e^{-{\omega_1 x^2\over 2}} e^{-{\omega_2 y^2\over 2}} $$
So there is no entanglement in the ground state between x and y. But if you look at it in a rotated basis (and $\omega_1 \ne \omega_2$), there is entanglement.
For a scalar quantum field in a spatial lattice in finite volume (time is still continuous), you have (if you Fourier transform in space) a bunch of decoupled harmonic oscillators (the sum on k is over non-redundant k's for a real scalar field, this is half the full space $k_x>0$):
$$ H = \sum_k {1\over 2} \dot{\phi_k}^2 + {k^2+m^2\over 2} \phi^2 $$
Which is a bunch of decoupled oscillators, so the ground state is;
$$ \psi_0(\phi_k) = \prod_k e^{-{\sqrt{k^2+m^2} |\phi_k|^2\over 2}} $$
That's not entangled in terms of $\phi_k$, but in terms of the $\phi_x$ (on the lattice), it is entangled. The vacuum wave-function Gaussian can be expressed here as:
$$ \psi_0(\phi) = e^{-\int_{x,y} \phi(x) J(x-y) \phi(y)} $$
Where $J(x-y) = {1\over 2} \sqrt{\nabla^2 + m^2} $ is not the propagator, it is this weird nonlocal square-root operator.
The vacuum for bosonic field theories is a statistical distribution, it is a probability distribution, which is the probability of finding a field configuration $\phi$ in a Monte Carlo simulation at any one imaginary time slice in a simulation (when you make the t-coordinate long). This is one interpretation of the fact that it is real and positive. The correlations in this probability distribution are the vacuum correlations, and for free fields they are simple to compute.
The axiomatic field theory material is not worth reading in my opinion. It is obfuscatory and betrays ignorance of the foundational ideas of the field, including Monte Carlo and path-integral.
General vacuum wave function for bosonic fields
In any path integral for bosonic fields with a real action (PT invariant theory), and this includes pure Yang-Mills theory and theories with fermions integrated out, the vacuum wave-function is the exact same thing as the probability distribution of the field values in the Euclidean time formulation of the theory. This is true outside of perturbation theory, and it makes it completely ridiculous that the rigorous mathematical theory doesn't exist. The reason is that the limits of probability distributions on fields as the lattice becomes fine are annoying to define in measure theory, since they become measures on distributions.
To see this, note that at t=0, neither the imaginary time nor the real time theory has any time evolution factors, so they are equivalent. So in an unbounded imaginary box in time, the expected values in the Euclidean theory at one time slice are equal to the equal time vacuum expectation values in the Lorentzian theories.
This gives you a Monte Carlo definition of the vacuum wave function of any PT invariant bosonic field theory, free or not. This is the major insight on ground states due to Feynman, described explicitly in the path integral and in the work on the ground state of liquid He4 in the 1950s (this is also a bosonic system, so the ground state is a probability distribution). It is used to describe the 2+1 Yang-Mills vacuum in the 1981 by Feynman (his last published paper), and this work is extended to compute the string tension by Karbali and Nair about a decade ago.
