Problem in derivation of nonlinear optics wave equation 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}$?
 A: You are correct, in general the term $\nabla \nabla\cdot \mathbf E$ should be kept in the equation.
The reason it is customary to drop it is that it simplifies the equation. The best justification I can think of is that this term should be negligible in isotropics if the susceptibility $\chi$ depends on the field only weakly.
Consider the case where the susceptibility does not vary at all, so $\chi(E)$ can be replaced by a constant independent of coordinates $\chi_0$. Since charge density is given by divergence of electric field and also by negative divergence of polarization, we have
$$
 \nabla \cdot \mathbf E = -\nabla \cdot (\chi_0 \mathbf E)
$$
and the only way this could hold everywhere in the isotropic medium$^*$ is if $\nabla \cdot \mathbf E = 0$ everywhere in the medium.
Purely transversal plane waves obey the condition $\nabla \cdot \mathbf E  =0.$
In those cases, the right-hand side of the wave equation does have two terms
$$
-\mu_0\frac{\partial^2\mathbf P}{\partial t^2} - \nabla \nabla\cdot \mathbf E
$$
but the second one is zero, so we can simplify the equation and work with RHS
$$
-\mu_0\frac{\partial^2\mathbf P}{\partial t^2}.
$$
Now, the effects of nonlinear optics are usually very small and thus dependence of $\chi$ on field or position is weak. Thus the second term $\nabla \nabla \cdot \mathbf E = 0$ , although non-zero, should be small with respect to the first term
$$-\mu_0\frac{\partial^2\mathbf P}{\partial t^2}.$$
We should still be able to neglect the second term and use the same equation as in the linear case. The nonlinearity will be preserved - but the variability of $\chi$ will enter the equation only through the first term, where polarization is nonlinear function of electric field strength.
In other words: since the terms in the equation are continuous functions of field and susceptibility, the small change in susceptibility should cause only small change in divergence of $\mathbf E$. The value should therefore be still negligible, at least in some range of field strengths that are not too high.
$*$ In crystals, the situation is different: because susceptibility is a tensor with differing components, zero divergence of $\mathbf D$ prevents zero divergence of $\mathbf E$ in all but special cases, such as when only one component of $\chi$ is realized (when polarization of the wave agrees with one of the principal axes of the crystal).
