Consider Maxwell's equations in a general medium without free charges or currents:
$$\nabla\cdot\vec{D} = 0 \tag{1}$$ $$\nabla\cdot\vec{B} = 0 \tag{2}$$ $$\nabla\times\vec{E} = -\frac{\partial \vec{B}}{\partial t} \tag{3}$$ $$\nabla\times\vec{H} = \frac{\partial \vec{D}}{\partial t} \tag{4}$$
The constitutive relations are $\vec{D} = \varepsilon_0\vec{E} + \vec{P}$ and $\vec{B} = \mu_0(\vec{H} + \vec{M})$.
Suppose for simplicity that the medium is nonmagnetic, so that $\vec{M} = \vec{0}$.
Additionally, suppose the medium is isotropic, homogeneous and nondispersive, so that we can write the nonlinear polarization density as $\vec{P} = \epsilon_0\chi(|\vec{E}|)\vec{E}$, where $\chi$ is (in general) a nonconstant scalar function of $|\vec{E}|$.
Taking the curl of equation $(3)$ and using $\vec{B} = \mu_0\vec{H}$ we get $$\nabla \times \nabla \times \vec{E} = -\mu_0\frac{\partial}{\partial t} (\nabla \times \vec{H}) = -\mu_0\frac{\partial^2 \vec{D}}{\partial t^2} \tag{5}$$
Now we can use the vector identity $\nabla \times \nabla \times \vec{E} = \nabla(\nabla \cdot \vec{E}) - \nabla^2\vec{E}$ and use the constitutive equation for $\vec{D}$ to obtain $$\nabla^2\vec{E} - \mu_0\varepsilon_0\frac{\partial^2 \vec{E}}{\partial t^2} = \mu_0\frac{\partial^2 \vec{P}}{\partial t^2} + \nabla(\nabla \cdot \vec{E}) \tag{6}$$
Now every source I have seen apparently assumes that $\nabla(\nabla\cdot \vec{E}) = \vec{0}$, after which they obtain the standard nonlinear optics wave equation with polarization source term: $$\nabla^2\vec{E} - \mu_0\varepsilon_0\frac{\partial^2 \vec{E}}{\partial t^2} = \mu_0\frac{\partial^2 \vec{P}}{\partial t^2} \tag{7}$$
Unfortunately, I am unable to see why this assumption is justified. The usual argument is that since $\vec{D} = \varepsilon_0(1 + \chi)\vec{E}$ then since $\nabla \cdot \vec{D} = 0$ we also have $\nabla\cdot\vec{E} = 0$.
But this argument is flawed in the nonlinear case since $\chi$ is in fact a function of $|\vec{E}|$, so that we need to use the product rule when computing the divergence operator: $$0 = \nabla\cdot\vec{D} = \varepsilon_0\nabla\chi(|\vec{E}|)\cdot\vec{E} + \varepsilon_0(1 + \chi(|\vec{E}|))\nabla\cdot\vec{E} \tag{8}$$
We therefore obtain $$\nabla\cdot\vec{E} = -\frac{\nabla\chi(|\vec{E}|)\cdot\vec{E}}{1 + \chi(|\vec{E}|)} \tag{9}$$
and it is far from obvious to me that this equals $0$ (or that the gradient of this expression equals $\vec{0}$) (indeed, $|\vec{E}|$ is generally spatially dependent so that the spatial derivatives are nonzero).
So is this assumption flawed and the most basic equation of nonlinear optics plainly wrong, or is there a correct way to derive the nonlinear optics wave equation (7), perhaps with additional assumptions on $\vec{E}$?