# Derivation of Hamiltonian in non-linear optics

Update 2: one half of the answer to this question can be found here and the other half in the second half of this post.

I asked this question a few hours ago, but decided to re-post it because I think I can phrase it more clearly. The question is closely related to this one and has to do with the derivation of a non-linear Hamiltonian for an electromagnetic field inside a dielectric, non-magnetic medium. It appears in a few places in the literature, for example: here, here, here, here and here.

It begins with the energy density of the EM field:

$$\begin{equation}\tag{1} \mathcal{H} = \int\boldsymbol{\mathrm{H}}\cdot d \boldsymbol{\mathrm{H}}+ \boldsymbol{\mathrm{E}}\cdot d \boldsymbol{\mathrm{D}}\end{equation}$$

Using $$\mu_0\boldsymbol{\mathrm{H}}=\boldsymbol{\mathrm{B}}$$ the first integrand is easy to evaluate, but the second is non-linear and a bit trickier. It can be written:

$$\begin{equation} \tag{2}\boldsymbol{\mathrm{E}}\cdot d \boldsymbol{\mathrm{D}} = \boldsymbol{\mathrm{E}} \cdot \bigg(\dfrac{\partial \boldsymbol{\mathrm{D}}}{\partial \boldsymbol{\mathrm{E}}} \bigg)d\boldsymbol{\mathrm{E}}\end{equation}$$

We know that: $$\begin{equation}\tag{3}\boldsymbol{\mathrm{D}}(\boldsymbol{\mathrm{E}})= \boldsymbol{\mathrm{P}}(\boldsymbol{\mathrm{E}})+\epsilon_0\boldsymbol{\mathrm{E}}\end{equation}$$

and can expand the polarisation density in a power series:

$$\begin{equation}\tag{4}\boldsymbol{\mathrm{P}} = \epsilon_0 (\chi^{(1)} \boldsymbol{\mathrm{E}}+ \chi^{(2)} \boldsymbol{\mathrm{E}}^{2} + \chi^{(3)} \boldsymbol{\mathrm{E}}^{3} + ... )\end{equation}$$

Now comes the step I'm having trouble justifying, the expression for the energy density obtained by integrating:

$$\begin{equation}\tag{5}\mathcal{H} = \epsilon_0\bigg(\frac{1+\chi^{(1)}}{2} \boldsymbol{\mathrm{E}}^{2} + \sum_{n\geq 2}^{N}\frac{n}{n+1}\chi^{(n)}\boldsymbol{\mathrm{E}}^{(n+1)} \bigg)+ \frac{1}{2\mu_0}\boldsymbol{\mathrm{B}}^{2}\end{equation}$$

specifically the scaling factors that appear for each term in the power series. Using the first term of $$\boldsymbol{\mathrm{E}}\cdot d \boldsymbol{\mathrm{D}}$$ as an example, we have:

$$\begin{equation}\tag{6}\boldsymbol{\mathrm{E}} \cdot \epsilon_0\dfrac{\partial \boldsymbol{\mathrm{E}}}{\partial \boldsymbol{\mathrm{E}}} d\boldsymbol{\mathrm{E}} = \epsilon_0\boldsymbol{\mathrm{E}} \cdot d\boldsymbol{\mathrm{E}}\end{equation}$$

and integrating this we get $$\frac{\epsilon_0}{2}\boldsymbol{\mathrm{E}}^2$$, which accounts for the first $$1/2$$ factor in the energy density. We also have the polarisation term:

$$\begin{equation}\tag{7}\dfrac{\partial \boldsymbol{\mathrm{P}}}{\partial \boldsymbol{\mathrm{E}}} d\boldsymbol{\mathrm{E}}\end{equation}$$

but looking at the first order terms, we have:

$$\begin{equation}\tag{8}\epsilon_0\bigg(\dfrac{\partial \boldsymbol{\mathrm{P}}^{(1)}}{\partial \boldsymbol{\mathrm{E}}}\bigg)_i = \chi^{(1)}_i\end{equation}$$

where $$i$$ represents both indices (getting a bug when I try to use two subscripts). Hence:

$$\begin{equation}\tag{9}\epsilon_0\boldsymbol{\mathrm{E}}\cdot \dfrac{\partial \boldsymbol{\mathrm{P}}^{(1)}}{\partial \boldsymbol{\mathrm{E}}} d\boldsymbol{\mathrm{E}} = \epsilon_0\boldsymbol{\mathrm{E}}\chi^{(1)}d\boldsymbol{\mathrm{E}} = \epsilon_0\sum E_i \chi_{i,k}^{(1)}dE_k \end{equation}$$

After integrating, the terms with $$i=k$$ will get a factor $$1/2$$, but the rest won't, leading to a discrepancy with the expression for $$\mathcal{H}$$.

If anyone is able to provide the correct reasoning, I'd be most grateful.

Update:

I believe I can partially answer my own question. First, I want to identify the mistake I made, which was in the integration. We have:

$$\begin{equation} \tag{10} \int \mathbf{E}\cdot \frac{\partial\mathbf{P}}{\partial\mathbf{E}} d\mathbf{E} \end{equation}$$

For the first order term: $$\begin{equation} \tag{11} \bigg( \frac{\partial\mathbf{P}^{(1)}}{\partial\mathbf{E}} \bigg)_{ij} = \chi^{(1)}_{ij} \end{equation}$$

and the integral is:

$$\begin{equation} \tag{12} \int E_i \chi^{(1)}_{ij} dE_j \end{equation}$$

The mistake was to evaluate the terms where $$i\neq j$$ as (setting $$i=x,j=y$$):

$$\begin{equation} \tag{13} \int E_x \chi^{(1)}_{xy} dE_y = E_x\chi_{xy}^{(1)}E_y \end{equation}$$

The above expression is wrong, because $$E_i$$ is not a variable, but a function of some implicit variable $$t$$, that defines the curve along which the integration of the field $$\mathbf{D}$$ is taken. The correct way of evaluating this is:

$$\begin{equation} \tag{14} \int E_x \chi^{(1)}_{xy} dE_y = \int_{t_0}^{t_1} E_x(t) \chi^{(1)}_{xy} \frac{dE_y(t)}{dt}dt = E_x(t)\chi_{xy}^{(1)}E_y(t)\bigg|_{t_0}^{t_1} - \int_{t_0}^{t_1} E_y(t) \chi^{(1)}_{xy} \frac{dE_x(t)}{dt}dt \end{equation}$$

If we take $$E_i(t_0)=0$$, $$E_i(t_1) = E_i$$ and use the fact that $$\chi^{(1)}$$ is symmetric under permutation of its indices, then: $$\begin{equation} \tag{15} \int E_x \chi^{(1)}_{xy} dE_y + \int E_y \chi^{(1)}_{yx} dE_x = E_x \chi^{(1)}_{xy} E_y \end{equation}$$

and this gives the correct first order contribution to the Hamiltonian density. We can also see it directly this way:

$$\begin{equation} \tag{16} \int E_i \chi^{(1)}_{ij} dE_j = \frac{1}{2} \bigg( \int E_i \chi^{(1)}_{ij} dE_j + \int E_j \chi^{(1)}_{ji} dE_i \bigg) = \frac{1}{2}\int d\big[\chi^{(1)}_{ij}E_iE_j\big] = \frac{1}{2}\chi^{(1)}_{ij }E_i E_j \end{equation}$$ where we used the symmetry again in the third step. For the second order term (I'm quite sure) that one has: $$\begin{equation} \tag{17} \bigg( \frac{\partial\mathbf{P}^{(2)}}{\partial\mathbf{E}} \bigg)_{ij} = 2\chi^{(2)}_{ijk}E_k \end{equation}$$ Using the same method as above, one finds: $$\begin{gather} \tag{18} 2\int E_i \chi^{(2)}_{ijk}E_j dE_k = \frac{2}{3} \bigg( \int E_i \chi^{(2)}_{ijk}E_j dE_k + \int E_i \chi^{(2)}_{ikj}E_k dE_j + \int E_j \chi^{(2)}_{jki}E_k dE_i \bigg) = \\ = \frac{2}{3} \int d\big[E_i \chi^{(2)}_{ijk}E_j E_k\big] = \frac{2}{3} E_i \chi^{(2)}_{ijk}E_j E_k \end{gather}$$ In the notation of the question, and cited papers, this is the same as: $$\begin{equation} \tag{19} \frac{2}{3}\chi^{(2)}\mathbf{E}^3 \end{equation}$$ It remains to justify the expression for the Jacobian in (17). For simplicity I considered a two-dimensional case, so the Jacobian would be: $$\begin{equation} \tag{20} \begin{bmatrix} a_{11} & a_{12}\\ a_{21} & a_{21} \end{bmatrix} \end{equation}$$ $$\begin{equation}a_{\gamma,\mu} = \frac{\partial P_{\gamma}^{(2)}}{\partial E_\mu} \end{equation}$$

and the first entry: $$\begin{gather} \tag{21} \frac{\partial P_{1}^{(2)}}{\partial E_1}= \frac{\partial}{\partial E_1} \big( \chi_{111}^{(2)} E_1 E_1 + \chi_{112}^{(2)} E_1 E_2 + \chi_{121}^{(2)} E_2 E_1 + \chi_{122}^{(2)} E_2 E_2 \big) = \\ = 2\chi_{111}^{(2)} E_1 + \chi_{112}^{(2)} E_2 + \chi_{121}^{(2)} E_2 = 2\big(\chi_{111}^{(2)} E_1 + \chi_{112}^{(2)} E_2\big) \end{gather}$$

where the last step is again justified by the symmetry of the tensor. Doing the same for all four entries one recovers the expression for the Jacobian in (17). I get a bit stuck trying to do this in index notation, because what I want to write is: $$\begin{equation} \tag{22} \bigg( \frac{\partial\mathbf{P}^{(2)}}{\partial\mathbf{E}} \bigg) = \sum_l \frac{\partial}{\partial E_l} \chi^{(2)}_{ijk} E_j E_k \end{equation}$$

that is, one should sum over all $$l$$, but this index shouldn't appear in the last three terms, so in normal index notation there would be no sum. This is also wrong: $$\begin{equation} \tag{23} \frac{\partial}{\partial E_j} \chi^{(2)}_{ijk} E_j E_k = 2\delta_{jk}\chi^{(2)}_{ijk}E_k + (1-\delta_{jk})\chi_{ijk}^{(2)}E_k \end{equation}$$ because the second expression doesn't get counted twice.

• Are you sure author's do not assume isotropic material? It would be strange to make such assumption since things like second harmonic generation are not allowed in isotropic materials, at the same time, without this assumption, your 5th equation (could you number them please) makes no sense. – Cryo Aug 25 at 23:05
• They do not assume that, but they do however assume that it is symmetric and real (because of assuming a lossless and non-magnetic material). I was thinking that maybe the susceptibility tensor is written in a diagonal form, since that would explain the factors in the sum (n from derivating, and 1/(n+1) from integrating again). So the question is if you can always diagonalise it under these assumptions. – fulis Aug 26 at 7:16
• Seems that the factors are an artifact of expression the Hamiltonian as integrals where the fields are the measures. Why do you need to do that? Why not express the energy density in terms of the fields and integrate it over space? That would give you the Hamiltonian with the unit of energy as expected. – flippiefanus Aug 26 at 7:33
• Your first reference it says " Here and below we keep implicit the usual tensor contraction of the susceptibilities with fields, the tensor χ (n) possessing n + 1 Cartesian components". I think this is where your Eq. (5) comes from. In which case, I think, when they say $\chi^{(n)}\mathbf{E}^{n+1}$, they mean $E^{\alpha'}g_{\alpha'\alpha}\chi^{(n)\alpha}_{\beta_1 \dots \beta_n}E^{\beta_1}\dots E^{\beta_n}$ – Cryo Aug 26 at 23:50
• Clearly these authors did not care about indices that much, or they would not try to skip over them. I would advise a better reference if I new one. For example, it makes no sense to me to work with even-order susceptibilities in this way. Either you do not have $\chi^{(2)}$, or you do have it, but then your space cannot be isotropic, but if your space is not isotropic you have to be very careful about all subsequent derivations. Afterall, it is a common trick in Physics to simply rotate the space in any way that simplifies your equations. – Cryo Aug 26 at 23:57