--------Definitions----------
For reference: Maxwell's equations as I understand them plus definitions of $\vec{B}$ and $\vec{D}$
\begin{align} \textrm{Maxwell's curl of }\vec{E} &: \nabla\times\vec{E} = -\partial_t\vec{B} \\ \textrm{Maxwell's curl of }\vec{H} &: \nabla\times\vec{H} = \partial_t\vec{E} + \vec{J} \\ \textrm{Gauss's law}&: \nabla\cdot\vec{D} = \rho \\ \textrm{No Magnetic Monopoles}&: \nabla\cdot\vec{B} = 0 \end{align}
Where assuming linear isotropic homogeneous media: \begin{align} \vec{D} = \epsilon\vec{E} \\ \vec{B} = \mu \vec{H} \end{align}
And $\rho$ is charge density and $\vec{J}$ is current density as usual.
--------Question-----------
When we convert the free-space Maxwell's equations into the Helmholtz equations, the procedure involves using the fact that $\nabla\cdot\vec{E} = \nabla \cdot \vec{B} = 0$. Specifically, it occurs when we take the curl-curl of the fields and simplify the results. So for the electric field for example we say:
\begin{align} \nabla \times \nabla \times \vec{E} &= \nabla(\nabla\cdot\vec{E}) - \nabla^2\vec{E} \\ &= -\nabla^2\vec{E} \end{align}
Where we've asserted $\nabla(\nabla\cdot\vec{E}) = 0$ by appealing to the divergenceless constraint on the electric field when we have no free charges.
Then, once we've used the above result with Maxwell's curl equation for the electric field, and done the same for the magnetic field using its own divergenceless constraint, and assumed a phasor form for the time dependence of the fields, we get to the Helmholtz equations which have simple plane wave solutions of constant velocity for all frequencies. For example, using Maxwell's curl equations for $\vec{E}$ and $\vec{B}$ in conjunction with the above identity gives:
\begin{align} \nabla\times\nabla\times\vec{E} = -\nabla^2\vec{E} &= -\partial_t\nabla\times\vec{B} \\ &=-\mu\epsilon\,\partial_t^2\vec{E} \end{align}
and assuming time dependence for $\vec{E}$ like $e^{-i\omega t}$ gives the Helmholtz equation for $\vec{E}$, which is just a 3D vector wave equation:
\begin{align} (\nabla^2 + \omega^2\mu\epsilon)\vec{E} = 0 \end{align}
A complete set of general solutions end up looking like plane waves (I know I'm being verbose here, sorry, but I want to be as clear as possible):
\begin{align} \vec{E}_{\vec{k},\omega}(\vec{r}, t) = \vec{E_0}e^{i(\vec{k}\cdot\vec{r} - \omega t)} \end{align}
with $\vec{k}$ and $\omega$ connected through the linear dispersion relation of velocity $c = \frac{1}{\sqrt{\mu\epsilon}}$ as usual.
But of course if we go back and check Gauss's law we see it's only satisfied for those $\vec{k}$ where $\vec{k}\cdot\vec{E}_0 = 0$. Thus, when we take the general solutions to the Helmholtz equations, we must reapply the divergenceless criteria of the fields to ensure that our fields are physical and satisfy Maxwell's equations in addition to the Helmholtz equation. It's this 'return' to to the divergenceless constraints (that we must satisfy after the fact by properly composing general solutions of the Helmholtz equations) that classically stipulates no longitudinal EM waves in free space.
But wait! Suddenly I'm confused. Didn't we already exhaustively apply all of Maxwell's equations in deriving the Helmholtz equation? Didn't we use the divergence equations to simplify our system to straightforward wave equations? How then can solutions pop out that defy the constraints we have implicitly deployed in the derivation of the very equations they solve? At what point did we 'loosen' the system of equations, so to speak - because it looks to me like we started with Maxwell and stayed with Maxwell the whole time? Why this need for postselection?
I suspect that this is just me conflating rules of solving equations algebraically with rules of solving differential equations. In some sense I understand that derivatives can erase information in an equation, so perhaps there's something (that to me is quite fuzzy) going on wen we take the curl of Maxwell's curl equations that erases somehow the divergence free constraints on the solutions that we nevertheless then immediately turn around and use to simplify the curl-curl identity.
I also understand that the divergence equations aren't independent of the curl equations, I believe so long as you axiomatically also assume the continuity equation and 'no magnetic monopoles'. If you take the divergence of Maxwell's curl equation for the magnetic field for example: \begin{align} \nabla\cdot(\nabla\times\vec{B}) = 0 = \partial_t\nabla\cdot\vec{E} + \nabla\cdot\vec{J} \end{align} Assuming the continuity equation - that $\nabla\cdot\vec{J} = -\partial_t \rho$, then gives Gauss's law for the electric field.
So perhaps therein lies the source of my confusion. We destroy the divergenceless condition when we take the curl of the curl equation? And even though we recycle it partially in simplifying the curl-curl identity, we have not sufficiently reconstrained the system?
I've been looking for a physics textbook or noteset (Griffiths, Jackson, Chew's notes on waveguides, ...) that explains this confusion of mine, but perhaps the issue is that it's the sort of 'senseless' confusion that doesn't warrant heading off. Maybe it's not much to do with physics and is really just a lack of understanding of PDE's and constraints on my part.
I'd be very grateful for any pointers in the right direction, whether to physics literature or math, that would settle this uneasiness I have. I hope I've explained myself well, please let me know if my confusion is confusing (I've tried to highlight the sentence I think most encapsulate my confusion).
Here is a question that abuts my own Are wave equations equivalent to Maxwell's equations in free space?. The top answer simply says 'no, you also need to use the divergence free constraints.' This is the sort of answer I keep seeing pop up, and I'm not able to logically placate the voice in my head that responds 'but didn't we already use them?'
