Non-linear optics Hamiltonian The non-linear optical Hamiltonian of the third order is written as:
$$
H = H_l - 
\frac{1}{3\varepsilon_0}\int d\mathbf{r} D^i(\mathbf{r})\Gamma_2^{ijk} D^j(\mathbf{r}) D^k(\mathbf{r}) - 
\frac{1}{4\varepsilon_0}\int d\mathbf{r} D^i(\mathbf{r})\Gamma_3^{ijkl} D^j(\mathbf{r}) D^k(\mathbf{r})D^l(\mathbf{r})
$$ where $H_l$ is the linear part. 
I have some problems getting this result; my approach is to take the  polarization 
$$
P_i = \Gamma^{ij}_1 D^j(\mathbf{r}) + \Gamma^{ijk}_2 D^j(\mathbf{r})D^k(\mathbf{r})+ \Gamma^{ijkl}_3 D^j(\mathbf{r})D^k(\mathbf{r})D^l(\mathbf{r})
$$ 
where the sum is implicit over repeated indexes and write the standard Hamiltonian as follow:
$$
H = \frac{1}{2}\int d\mathbf{r} \mathbf{E}\cdot\mathbf{D} = \frac{1}{2\varepsilon_0}\int d\mathbf{r} (\mathbf{D-P})\cdot\mathbf{D} = \frac{1}{2\varepsilon_0}\int d\mathbf{r}\mathbf{D}\cdot\mathbf{D} - \frac{1}{2\varepsilon_0}\int d\mathbf{r}\mathbf{P}\cdot\mathbf{D}\, .
$$ 
(Taken from section VII of this paper.)
Focusing on the last term:
$$\int d\mathbf{r}\mathbf{P}\cdot\mathbf{D} = \int d\mathbf{r}P^iD^i =\int d\mathbf{r}D^i(\mathbf{r})\left(\Gamma^{ij}_1 D^j(\mathbf{r}) + \Gamma^{ijk}_2 D^j(\mathbf{r})D^k(\mathbf{r})+ \Gamma^{ijkl}_3 D^j(\mathbf{r})D^k(\mathbf{r})D^l(\mathbf{r})\right).
$$ 
This leads to a very similar Hamiltonian apart from the factors $\frac{1}{3},\frac{1}{4}$ in front of the integral. What am I doing wrong? Is there an explanation for those factors?
 A: Your expression for the hamiltonian, $H = \frac12\int\mathbf E\cdot \mathbf D \, \mathrm d\mathbf r$, is incorrect. Following up the reference trail of the paper you quote quickly leads to a related paper that makes it clear that the relevant quantity in that formalism is
$$
H = \int\left[ \int \mathbf E\cdot \mathrm d\mathbf D\right]\mathrm d\mathbf r,
$$
and further pointing to chapter 4 of J.D. Jackson's Classical Electrodynamics as a reference for that. 
In principle, the factors of $1/3$ and $1/4$ come out of the integration with respect to $\mathbf D$. Thus, for the linear case, the relevant integral is of the form
$$
\int D^j\mathrm dD^i=\frac{2-\delta^{ij}}{2}D^jD^i,
$$
with a prefactor to handle the behaviour where $i=j$ and therefore the integral should read e.g. $\int D^1\mathrm dD^1=\frac{1}{2}(D^1)^2$. As you can imagine, the higher-order terms get messier on  a combinatorically fast scale, because of the many ways the polynomial could integrate; as an example, the second-order nonlinearity gives a term of the form
$$
\int D^jD^k\mathrm dD^i=\left(1-\frac{\delta^{ji}+\delta^{ki}}{2}+\frac{\delta^{ji}\delta^{ki}}{3}\right)D^jD^kD^i 
$$
to account for the cases where $j=i$ or $k=i$ separately, and where $i=j=k$. 
Finally, you need to contract this with the $\Gamma_2^{ijk}$, using the fact that they are symmetric in the Cartesian components (because the line integral $\int \mathbf E\cdot \mathrm d\mathbf D$ needs to be independent of the path of integration). This should give you the results that Sipe et al. claim in that second paper, but to be honest with you I'm not sure this actually works, even in the linear case: the integral still gives
$$
\int \Gamma_1^{ij}D^j\mathrm dD^i=\frac{2-\delta^{ij}}{2}\Gamma_1^{ij}D^jD^i=\frac12\Gamma_1^{xx}D^xD^x+(\Gamma_1^{xy}+\Gamma_1^{yx})D^xD^y+\frac12\Gamma_1^{yy}D^yD^y
$$
when only the $x$ and $y$ components are nonzero, and this does not reduce to $\Gamma_1^{ij}D^jD^i$ as claimed. Ultimately this can be fixed by refactoring those combinatoric terms into the $\Gamma$s, but they seem to be explicitly not doing that, so I would tread with extreme caution with those results.
