What is a vectorial vortex beam? I know that a Laguerre-Gaussian (LG) beam with topological phase non zero has helical wavefront, but still linear polarization. So if LG beams are not vectorial vortices, what could they be?
 A: Generally "vector" and "scalar" are models for the same kind of thing, the former can be more accurate than the latter. Generally, "vector" models of the electromagnetic field tend to be needed for "fast" fields: i.e. those of high numerical aperture and which contain a wide angular spread of wave directions. "Vector" models are also needed for fields whose polarisation significantly varies across the transverse plane. Electromagnetism is inherently a vector phenomenon (i.e. has polarisation), so the words "vector" and "scalar" don't reflect fundamentally different fields, they refer to which approximation one can get away with. Even the linear polarised LG beam cannot have a perfectly uniform polarisation - the only field that can have this is a plane wave - so, depending on the accuracy needed in simulation, you may have to model this field as a fully vector phenomenon. The "faster" the LG field is, the more inaccurate is the scalar model. So an LG beam of the kind you describe with a Gaussian waist parameter of, say, half a wavelength will be significantly "vector" and you'll see significant inaccuracies when modelling it as a scalar field.
In a homogeneous medium, each Cartesian component of the EM field (and of the Lorenz gauged 4-vector) fulfills the Helmholtz equation. So, if the polarisation is mostly in one transverse direction, one scalar field $\psi = E_x(\vec{r},\,t)$ can define the whole EM  field. That is, $E_y\approx 0$, whence you can, with the assumption of unidirectional propagation (i.e. all field components propagating with the same sign $z$ (axial) component) retrieve the $E_z$ component from Gauss's law $\nabla \cdot \vec{E}=0$ (i.e. the field is solenoidal). In Fourier (momentum) space, Gauss's law becomes the statement that the transverse components are orthogonal to the wavevector, so if there is no $y$ component, the $\tilde{E}_x(\vec{k})$ in Fourier space wholly defines the $z$ component through $k_x\,\tilde{E}_x(\vec{k}) + k_z \,\tilde{E}_z(\vec{k})=0$ and $k_z=\sqrt{k^2-k_x^2}$ (here's where you need unidirectionality: so you know the sign of the square root). Once you have the whole $\vec{E}$ field, you can calculate $\vec{H}$ through Faraday's law and you're done. Thus the field is accurately described by the lone scalar. In the more complicated case of near uniform circular polarisation, you might think this becomes needfully "vector", but not so. If there are mostly left or mostly right hand components, you can work the same trick with the Riemann-Silberstein vectors see my answer here for more information, so the problem stays essentially "scalar".
To overcome this problem, one in general needs to model and track two field components: reasoning as above, $E_x$ and $E_y$ will be enough to define the whole, unidirectionally propagating field. The circular polarised analogy is that you track the components $F_{\pm} = E_x \pm\,i\,c\,B_y$, which yields two scalars standing for the left and right handed amplitudes. This is the approach used, at least to some degree, in the "LightTrans Virtual Lab" optical modelling software, although I have no experience of this beast.
