Scalar wave approximation to electrodynamic wave equation
The art of going from scalar optics to a full solution of the Maxwell equations is known as vector optics; it is a very wide field with multiple different methods, and ultimately which one you choose depends on the situation. Ultimately, the underlying theme is your
Question: Given a solution $u(\mathbf r)$ of the scalar wave equation $\frac{\partial^2u}{\partial t^2}=c^2\nabla^2 u$, how can one obtain a full set of vector solutions to the Maxwell equations?
As you note, taking $\mathbf E(\mathbf r) = \mathbf r \, u(\mathbf r)$ is not going to cut it, and neither is the constant-polarization solution $\mathbf E(\mathbf r) = \mathbf E_0 u(\mathbf r)$, which will fail to have a zero divergence if $u(\mathbf r)$ has any dependence along $\mathbf E_0$, so you need to look elsewhere.
For this, there's a number of things you can try. (Most of these assume that you start with a solution of the Helmholtz equation $(\nabla^2 +k^2)u=0$ rather than with full time dependence, but you can roll back from there through a Fourier transform if you really need to.) So, some techniques to manufacture vector solutions of the Helmholtz equation out of known solutions:
- Multiply a scalar solution $u(\mathbf r)$ by a constant polarization $\mathbf E_0$ to get $\mathbf E(\mathbf r) = \mathbf E_0 u(\mathbf r)$.
- Take the gradient $\mathbf E(\mathbf r) = \nabla u(\mathbf r)$ of a scalar solution $u(\mathbf r)$.
- Take the cross product of an existing vector solution $\mathbf E_0(\mathbf r)$ with some constant vector $\mathbf c$ to get $\mathbf E(\mathbf r)= \mathbf c\times\mathbf E_0(\mathbf r)$.
- Take the curl of an existing vector solution $\mathbf E_0(\mathbf r)$ to get $\mathbf E(\mathbf r)= \nabla \times\mathbf E_0(\mathbf r)$.
Depending on the situation, most of these won't work, but some of these might work (and indeed, you might need to daisy-chain several of these methods to get the solutions you want; thus, one set of lecture notes in my back archive uses things like $\nabla\times\nabla\times\left(\mathbf c \, u(\mathbf r)\right)$ as one of its solutions). Choosing a method that will give you a solution that is simple to work with and appropriate to your problem can be a bit of an art, and again it depends on the situation. Similarly, depending on the regime (like, say, a loosely focused gaussian beam), several of the techniques above which look like they don't work (like a constant polarization) can work rather well as approximations - but only on a case-by-case basis.
What this means, ultimately, is that the relationship between the solutions of vector optics and the simpler solutions of scalar optics will depend on the problem, and there is no one-size-fits-all description of that relationship. This then percolates to the fact that generally speaking $u(\mathbf r)$ is not an approximation for $\mathbf E(\mathbf r)$: instead, it is an auxiliary quantity on the way up the ladder towards a full vector solution of the Maxwell equations.
To be a bit more precise, $u(\mathbf r)$ isn't just an auxiliary quantity ─ it's usually something that you can derive from the full vector solution, normally by just taking one of its components, since if you have a full vector solution $\mathbf E(\mathbf r)$ then all of its cartesian components will be wave-equation solutions. That means that while $u(\mathbf r)$ isn't sufficient to construct a $\mathbf E(\mathbf r)$, it is in a sense necessary, and often we can get much of the relevant insight about our desired vector solutions just from the behaviour of the scalar solutions in the geometry of interest. Thus, even when you know that it is only a partial description of your system, it is still useful to study the dynamics of scalar waves in whatever geometry you're working.
On the other hand, however well you understand the scalar-wave dynamics, you still need to construct a full vector solution, and if you want a rough-and-ready most-common-case answer for how this gets done, though, it turns out to be the first one above, which you mentioned in the question: the constant-polarization solution $$\mathbf E(\mathbf r) = \mathbf E_0 u(\mathbf r).$$ It is easy to see that this is a solution to the wave equation, but as you say,
this could violate badly the divergence free condition.
However, "badly" is a relative term, and - particularly if the characteristic dimension $a$ of the problem are large compared to the wavelength - it can work just fine as an approximation. If you really want a 'true' solution in those cases, you can find them out as above, and then you recover the scalar-optics constant-polarization solution as the limit $\lambda\ll a$ of the true solutions. (In this connection, see my answer to this question and its references for examples of vector optics in action.) However, most often, it's really not necessary to go through all that hassle, because the scalar-optics solution works perfectly fine to explain experiment, and that is the ultimate goal.