Let's say a gaussian beam propagates along z-axis. According to Wikipedia its intensity will be: $I(r,z) = I_0 \frac{\omega_0}{\omega(z)}\exp(\frac{-2r^2}{\omega(z)^2})$
I think of the intensity as a Poynting vector. In the definition above it is a scalar.
Question: is it assumed that in the paraxial approximation, the energy flows predominantly in the z-direction? That is to say, Poynting vector is parallel to the z-axis.
You should not think about intensity as a Poynting vector because light intensity, usually measured in W/m$^2$, is a scalar value. However, the time-averaged magnitude of a Poynting vector is the intensity, so both are indeed related.
The paraxial approximation is not assuming energy flow only in the $z$ direction, after all, it is giving us inherently divergent Gaussian beams. One can even argue, that the presented formula is exact if the energy flux is related to the surface normal to your Poynting vector at a given point. On the other hand, if the intensity is defined on a flat surface, then there is an assumption that your beam consists of planar waves.
