The passage from the full time-dependent wave equation $(\mathrm{W})$ to the Helmholtz equation $(\mathrm{H})$ is nothing more, and nothing less, than a Fourier transform. For sufficiently regular functions, both $u$ and $F$ can be written as superpositions of monochromatic fields, i.e. they are both Fourier transforms of the form
$$
u(x,t) = \int_{-\infty}^\infty U(x,\omega) e^{-i\omega t} \mathrm d\omega
\quad\text{and} \quad
f(x,t) = \int_{-\infty}^\infty F(x,\omega) e^{-i\omega t} \mathrm d\omega,
$$
where the temporal Fourier coefficients $U(x,\omega)$ and $F(x,\omega)$ depend on the position - or, switching perspectives, they give us functions of $x$ for each $\omega$. Now, all we've done so far is a fancy rewriting of our variables, but there are two crucial aspects of the wave equation that make this useful:
- it is linear, and
- the only dependence on time is through $\partial_t^2$, which is a linear operator whose eigenfunctions are precisely the Fourier kernel, i.e. $\partial_t^2 e^{-i\omega t} = -\omega^2 e^{-i\omega t}$.
The linearity allows us to break in the wave equation's linear operators all the way through to the Fourier coefficients, and the eigenvalue relation for $\partial_t$ enables us to switch that partial differentiation to an algebraic factor on that sector, giving us
\begin{align}
0
& =
-\partial_{t}^2 u(x,t) + c^2 \nabla^2 u(x,t) + f(x,t)
\\ & =
-\partial_{t}^2 \int_{-\infty}^\infty U(x,\omega) e^{-i\omega t} \mathrm d\omega
+ c^2 \nabla^2 \int_{-\infty}^\infty U(x,\omega) e^{-i\omega t} \mathrm d\omega
+ \int_{-\infty}^\infty F(x,\omega) e^{-i\omega t} \mathrm d\omega
\\ & =
\int_{-\infty}^\infty \left[
-U(x,\omega) \partial_{t}^2 e^{-i\omega t}
+ c^2 \nabla^2 U(x,\omega) e^{-i\omega t}
+ F(x,\omega) e^{-i\omega t}
\vphantom{\sum}\right]\mathrm d\omega
\\ & =
\int_{-\infty}^\infty \left[
\omega^2U(x,\omega)
+ c^2 \nabla^2 U(x,\omega)
+ F(x,\omega)
\vphantom{\sum}\right] e^{-i\omega t} \mathrm d\omega
. \tag{1}
\end{align}
From here, it's easy to see that if $f(x,t)$ is given (so $F(x,\omega)$ is also given), we can find a solution of the original equation by setting $U(x,\omega)$ to be a solution of the Helmholtz equation,
$$
c^2 \nabla^2 U(x,\omega) + \omega^2U(x,\omega) = - F(x,\omega)
$$
or with the cosmetic change $k=\omega/c$,
$$
\nabla^2 U(x,\omega) + k^2U(x,\omega) = - \frac{1}{c^2} F(x,\omega).
$$
(In addition, it's easy to show that the Fourier transform in $(1)$ means that this is a necessary condition, but if all you're doing is finding solutions, as opposed to characterizing the general solution, then the sufficiency is enough.)
OK, so that is the formal side. How is this used in the real world? There are three main ways that one uses this.
The clearest is when the wave equation is being forced by a source that is itself monochromatic (or close enough to monochromatic that your situation doesn't care about the difference), or in terms of the Fourier amplitude $F(x,\omega) = F(x) \delta(\omega-\omega_0)$. This is the case, for example, when one considers the electromagnetic emission of an antenna set to a very narrow band of frequencies.
Physically speaking, the Helmholtz equation $(\mathrm{H})$ does encode propagation, in a very real sense ─ except that you're considering in one single go the coherent superposition of the emission that comes from a source that is always turned on, and oscillating at a constant frequency for all time. In this case, you expect the physical response to be at that same frequency, but the spatial response can be complicated in the presence of reflections, dispersive media, or whatnot; we solve the Helmholtz equation to find that spatial response.
In this case, $\omega$ is obviously fixed by the external driver.
A separate application is when we solve for resonant modes of the domain in question; these are nonzero solutions to the Helmholtz equation that hold even when the driver $F$ is zero, and they are important e.g. when $F$ is an impulse that's confined in time, like hitting a drum, and the effects are left to resonate in a confined domain which the energy cannot leave easily.
In this case, $\omega$ is fixed by the domain to one of a discrete set of resonant frequencies that sustain nonzero solutions of $(\mathrm{H})$ even when the driver is zero, and the process of solving the Helmholtz equation includes finding those resonant frequencies. The end goal in this calculation is a set of resonant frequencies $\{\omega_n\}$ with a corresponding set of solutions $\{u_n(x)\}$ which satisfy the homogeneous Helmholtz equation at that frequency and which form a complete basis, in the $L_2$ sense, for functions over the domain in question.
To reconcile this with the driver, the simplest case is to consider an impulsive driver, i.e. something of the form $f(x,t) = f(x)\delta(t)$, with a flat Fourier transform. In this case, all modes see the impulse, but only the resonant modes are able to respond. In this case, you decompose $f(x)$ as a linear combination of the $u_n(x)$, and this tells you how much each mode gets excited, which determines the temporal evolution after the impulse is gone.
Does this describe "propagation" in a suitable sense? Well, you're ultimately solving for the propagation of an initial impulsive disturbance, like plucking a string, by finding a clever decomposition of that initial disturbance in terms of modes that evolve cleanly (monochromatically) in time. So, yes.
Finally, there is also the case where you just have some arbitrary driver $f(x,t)$ for the wave equation, and all that you can say about its Fourier transform is that it exists. In this case, $\omega$ isn't 'chosen', as such: instead, it is a continuous parameter of the problem, where you solve a continuous set of separate inhomogeneous Helmholtz equations to get the $U(x,\omega)$, and then you add them all up coherently to get $u(x,t) = \int_{-\infty}^\infty U(x,\omega) e^{-i\omega t} \mathrm d\omega$.
It's important, however, not to underestimate the important of what you can say about $F$: just by saying "the temporal Fourier transform of $f(x,t)$ exists", you're saying that $f(x,t)$ can be understood as a superposition of monochromatic waves, each of which can be solved independently and which will cause some monochromatic response $U(x,\omega) e^{-i\omega t}$, which can then be added together to give the global response to the driver.
Best Answer
Yes, indeed you can use your knowledge of the scalar Helmholtz equation. The difficulty with the vectorial Helmholtz equation is that the basis vectors $\mathbf{e}_i$ also vary from point to point in any other coordinate system other than the cartesian one, so when you act $\nabla^2$ on $\mathbf{u}$ the basis vectors also get differentiated. This forces you to calculate $\nabla^2 \mathbf{u}$ through the identity $$ \nabla^2 \mathbf{u} = \boldsymbol{\nabla} (\boldsymbol{\nabla} \cdot \mathbf{u}) - \boldsymbol{\nabla}\times (\boldsymbol{\nabla}\times \mathbf{u}) \tag{1} $$ which is really cumbersome to deal with by brute force. A smart way to avoid all the hassle is by using the ansatz $$ \mathbf{u} = \mathbf{r} \times (\boldsymbol{\nabla} \psi) \tag{2} $$ where $\psi$ satisfies the scalar Helmholtz equation $$ (\nabla^2 + k^2) \psi = 0. $$
To check that $(\nabla^2 + k^2) \mathbf{u} = 0$ yourself you have to plug the ansatz $(2)$ on $(1)$ and make use of many vector identities and the scalar Helmholtz equation. The calculation is quite involved, so I'll point you to check Reitz, Milford & Christy's Foundations of Electromagnetic Theory, there they do the full calculation. With ansatz $(2)$ proven, it's just a matter of plugging the relevant mode $\psi_{lm}$ in eq. $(2)$ that you get your solution $\mathbf{u}_{lm}$.