Short answer: $\nabla \cdot \mathbf{D}=0$ comes from taking the divergence of $\nabla \times \mathbf{H}_F = j \omega \epsilon \mathbf{E}_F$ since we have assumed that $\mathbf{J}=0$ for this case.
We use surface equivalence to replace physical structures with fictitious $\mathbf{J}$ and $\mathbf{M}$ so that we can use the free space Green's function to calculate radiated electric and magnetic fields.
And now to overexplain:
For the case of both electric and magnetic sources, we start with the equations
$\nabla \times \mathbf{E} = -j\omega \mu \mathbf{H} -\mathbf{M}$
$\nabla \times \mathbf{H} = j\omega \epsilon \mathbf{E} + \mathbf{J}$
Using superposition, we can consider the electric and magnetic current sources separately.
$\nabla \times \mathbf{E}_A = -j\omega \mu \mathbf{H}_A$
$\nabla \times \mathbf{H}_A = j\omega \epsilon \mathbf{E}_A + \mathbf{J}$
and
$\nabla \times \mathbf{E}_F = -j\omega \mu \mathbf{H}_F - \mathbf{M}$
$\nabla \times \mathbf{H}_F = j\omega \epsilon \mathbf{E}_F$
and using superposition we have
$\mathbf{E} = \mathbf{E}_A+\mathbf{E}_F$
$\mathbf{H} = \mathbf{H}_A+\mathbf{H}_F$
Solving the electric current case, we use the usual magnetic vector potential $\mathbf{A}$, where $\mathbf{H}_A = \frac{1}{\mu}\nabla \times \mathbf{A}$, and this will yield $\mathbf{E}_A$ and $\mathbf{H}_A$. All of the logic for solving for $\mathbf{E}_F$ and $\mathbf{H}_F$ from the magnetic current source $\mathbf{M}$ is exactly the same, just switch $\mathbf{E}$ and $\mathbf{H}$ (and other constants and signs) and do the same thing.
Take the following divergence
$\nabla \cdot \left( \nabla \times \mathbf{H}_F = j\omega \epsilon \mathbf{E}_F \right)$
to get
$\nabla \cdot \epsilon \mathbf{E}_F=0$
I think this was the part you were missing. Then we can write
$\mathbf{E}_F = -\frac{1}{\epsilon} \nabla \times \mathbf{F}$.
This will allow us to solve for $\mathbf{E}_F$ and $\mathbf{H}_F$.
Leaving some steps out, we'll end up with the vector potentials $\mathbf{A}$ and $\mathbf{F}$ as
$\nabla^2 \mathbf{A} +k^2 \mathbf{A} = -\mu \mathbf{J}$
$\nabla^2 \mathbf{F} +k^2 \mathbf{F} = -\epsilon \mathbf{M}$
This can be seen as 6 separate scalar wave equations, and in free space the solution is well known using Green's functions. Note that using the Lorenz gauge is necessary to obtain this form. Assuming that $\mathbf{J}$ and $\mathbf{M}$ are surface currents on some surface $\Gamma$, then
$\mathbf{A}(\mathbf{r}) = \mu \int_\Gamma g(\mathbf{r},\mathbf{r}') \mathbf{J}(\mathbf{r}') ds'$
$\mathbf{F}(\mathbf{r}) = \epsilon \int_\Gamma g(\mathbf{r},\mathbf{r}') \mathbf{M}(\mathbf{r}') ds'$
where $g(\mathbf{r},\mathbf{r}') = \frac{e^{-jk|\mathbf{r}-\mathbf{r}'|}}{4\pi|\mathbf{r}-\mathbf{r}'|}$
And then once we have $\mathbf{A}(\mathbf{r})$ and $\mathbf{F}(\mathbf{r})$ we can get
$\mathbf{E}(\mathbf{r}) = -j \omega \mathbf{A} - \frac{j}{\omega \mu \epsilon} \nabla (\nabla \cdot \mathbf{A}) - \frac{1}{\epsilon} \nabla \times \mathbf{F}$
$\mathbf{H}(\mathbf{r}) = -j \omega \mathbf{F} - \frac{j}{\omega \mu \epsilon} \nabla (\nabla \cdot \mathbf{F}) + \frac{1}{\mu} \nabla \times \mathbf{A}$
We can also take the limit $|\mathbf{r}| \rightarrow \infty$ to obtain far fields.
As to why we consider fictitious magnetic currents $\mathbf{M}$:
Suppose we have a radiating antenna, and we know the electric and magnetic fields in some close vicinity of the antenna (e.g. using a numerical finite element simulation), but we want to find the radiated fields at some far away point. We can't use simply use the actual electric currents along with the integrals above, since those equations rely on the free space Green's function and the presence of a physical radiating structure means this is not free space.
Given a closed surface $\Gamma$ which contains the antenna, if we know the fields $\mathbf{E}$ and $\mathbf{H}$ on $\Gamma$, then we can consider a equivalent problem in which we (1) introduce fictitious sources $\mathbf{J}=\mathbf{\hat{n}} \times \mathbf{H}$, $\mathbf{M} = -\mathbf{\hat{n}} \times \mathbf{E}$ on $\Gamma$, where $\mathbf{E}$ and $\mathbf{H}$ are the known fields from the original problem, and (2) remove the physical structure and set $\mathbf{E}=\mathbf{H}=0$ inside the surface. Surface equivalence tells us that outside $\Gamma$ the two scenarios will produce the same fields. The equivalent problem has the advantage that it is free space, and therefore we know the Green's function. Therefore we can use the integrals above to calculate the vector potentials and then the electric and magnetic fields at any point outside $\Gamma$.