## Classical origin of optical nonlinearity The classical origin of linear susceptibility is from [[Classical Lorentz model]], namely model the electron behaviors as forced damping. The bounding force is considered as spring forces, namely parabolic potentials. We could, using the same manner, explain the origin of higher order susceptibility. To put it at the beginning, the origin of 2nd (or any even) order nonlinearity is by non-centrosymmetric potentials. For odd order nonlinearity, centrosymmetric potential (but non-parabolic) acts a origin. One shortcoming of such classical model is there is **only one resonance frequency $\omega_{0}$**. Therefore, the full picture still requires quantum mechanical description, who processes more than one energy eigenvalues. But the classical model still provides good estimations **when optical frequencies (namely input frequencies) are considerably smaller than the lowest electronic resonance frequency of the material system**. >[!Notice]- >- A medium is **centrosymmetric** if its crystal structure or molecular arrangement has a **center of symmetry**. This means that for every point $(x, y, z)$ in the material, there is an equivalent point $(-x, -y, -z)$ where the properties of the medium are identical. >- A medium is **non-centrosymmetric** if its crystal structure or molecular arrangement lacks a center of symmetry. This absence of inversion symmetry allows for both even-order and odd-order nonlinear optical processes. >- It is the symmetry properties of atoms/molecules make the electron potential energy function have the corresponding symmetry. ### Non-centrosymmetric media #### System setup Like [[Classical Lorentz model|Lorentz model]], we description the force of an electron feels and write the Newton's second law. Because we did not write the symmetry properties of the material and how the material is placed in the field, we could not show the tensor nature of the $\chi$. For simplification, we use scalar to write these quantities. We have $\colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 \tilde{x}}{\mathrm{d}t^2}\)}+\colorbox{#ffa94d66}{\(2 \gamma \frac{\mathrm{d}\tilde{x}}{\mathrm{d}t}\)}+\colorbox{#4dff5062}{\(\omega_0^2 \tilde{x} + a \tilde{x}^2\)}=\colorbox{#ff89a962}{\(-\frac{e \tilde{E}(t)}{m}\)}$ The $\colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 \tilde{x}}{\mathrm{d}t^2}\)}$ part is the acceleration. The second part $\colorbox{#ffa94d66}{\(2 \gamma \frac{\mathrm{d}\tilde{x}}{\mathrm{d}t}\)}$ is the damping force, $F_{damping}/m$. It is worth to note that here we have the damping coefficient $2\gamma$, which makes the this and later expression for $D$, the complex denominator function looks different from [[Classical Lorentz model#Polarization and relative permittivity under Lorentz model]]. The third part $\colorbox{#4dff5062}{\(\omega_0^2 \tilde{x} + a \tilde{x}^2\)}$ is from the electronic potential, in original model writes as $F_{spring}/m$. But here we have the $a \tilde{x}^2$ term, where $a$ is a parameter that characterize the strength of nonlinearity, and we call it the restoring force. And finally, $\colorbox{#ff89a962}{\(-\frac{e \tilde{E}(t)}{m}\)}$, is the $F_{driving}/m$. the electric field. The `tilde` still means oscillation with time. Here we did not the proper sign of these quantities. Write the restoring force, namely the $\colorbox{#4dff5062}{green}$ part, ${\tilde{F}}_{\text{restoring}} = -m \omega_0^2 \tilde{x} - ma \tilde{x}^2$ and the potential would be, $U(\tilde{x}) = - \int \tilde{F}_{\text{restoring}} \, d\tilde{x} = \frac{1}{2} m \omega_0^2 \tilde{x}^2 + \frac{1}{3} ma \tilde{x}^3$ Namely, a parabolic potential plus a non-parabolic potential. This is non-centrosymmetric since $U(\tilde{x})\neq U(-\tilde{x})$. And it looks like the following plot, ![[Drawing 2024-09-09 15.25.31.excalidraw.svg]] Now with the damping model set, we may define our optical (electric) field, $\tilde{E}(t) = E_1 e^{-i \omega_1 t} + E_2 e^{-i \omega_2 t} + \text{c.c.}$ where $E_{1}=E(\omega_{1})$, $E_{2}=E(\omega_{2})$. Since we mainly concern 2nd order nonlinearity, two frequencies should be enough. Although we remove the time dependence, these $E$ are still complex functions containing position related phase. Now the problem is to solve ${\frac{\mathrm{d}^2 \tilde{x}}{\mathrm{d}t^2}}+{2 \gamma \frac{\mathrm{d}\tilde{x}}{\mathrm{d}t}}+{\omega_0^2 \tilde{x} + a \tilde{x}^2}={-\frac{e \tilde{E}(t)}{m}}$ #### Solve the problem with perturbation ##### First order Define a perturbation parameter $\lambda$ and multiple it to the field, ranges from $0$ to $1$ and be one at the end of the calculation. The the problem becomes ${\frac{\mathrm{d}^2 \tilde{x}}{\mathrm{d}t^2}}+{2 \gamma \frac{\mathrm{d}\tilde{x}}{\mathrm{d}t}}+{\omega_0^2 \tilde{x} + a \tilde{x}^2}={-\lambda\frac{e \tilde{E}(t)}{m}}$ Expand $x$ in a power series form, we have $\tilde{x} = \lambda \tilde{x}^{(1)} + \lambda^2 \tilde{x}^{(2)} + \lambda^3 \tilde{x}^{(3)} + \cdots$ Then reorganize with respect to the order of $\lambda$, namely $\lambda$, $\lambda^{2}$, $\lambda^{3}$, $\cdots$, we have $\lambda$: $\frac{\mathrm{d}^2 \tilde{x}^{(1)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(1)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(1)} = -\frac{e \tilde{E}(t)}{m}$ $\lambda^{2}$: $\frac{\mathrm{d}^2 \tilde{x}^{(2)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(2)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(2)} + a \left[ \tilde{x}^{(1)} \right]^2 = 0$ $\lambda^{3}$: $\frac{\mathrm{d}^2 \tilde{x}^{(3)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(3)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(3)} + 2a \tilde{x}^{(1)} \tilde{x}^{(2)} = 0$ $\cdots$ The equation with $\lambda$ is exactly the classical linear Lorentz model. We can directly give the result, $\tilde{x}^{(1)}(t) = x^{(1)}(\omega_1) e^{-i \omega_1 t} + x^{(1)}(\omega_2) e^{-i \omega_2 t} + \text{c.c.}$ $x^{(1)}(\omega_j) = -\frac{e}{m} \frac{E_j}{\omega_0^2 - \omega_j^2 - 2i \omega_j \gamma}$ Or write as the simplified form with the complex denominator function, $x^{(1)}(\omega_j) = -\frac{e}{m} \frac{E_j}{D(\omega_j)}$ $D(\omega_j) = \omega_0^2 - \omega_j^2 - 2i \omega_j \gamma$ As stated before, the $2\gamma$ is from our setup, we have $2\gamma$ as the damping coefficient. This gives the overall displacement as a superposition of each input frequency. ##### Second order And to gain the second order response, one should check the $\lambda^{2}$ expression, who contains $[x^{(1)}]^{2}$. This will automatically provide $2\omega$ term by $[e^{i\omega t}]^{2}$. For example, if we wanna get the response at $2\omega_{1}$ substitute the $x^{(1)}(\omega_{1})e^{-i\omega_{1} t}$ and $x^{(1)}(\omega_{1})e^{-i\omega_{1} t}$, we have $\frac{\mathrm{d}^2 \tilde{x}^{(2)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(2)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(2)} = \frac{-a \left( \frac{e E_1}{m} \right)^2 e^{-2i \omega_1 t}}{D^2(\omega_1)}$ And we are looking for a solution with the form $\tilde{x}^{(2)}= x^{(2)}(2\omega_{1})e^{-2i \omega_{1}t}$, substitute this into the above equation, we have $x^{(2)}(2\omega_1) = \frac{-a \left( \frac{e}{m} \right)^2 E_1^2}{D(2\omega_1) D^2(\omega_1)}$ The $D$ has the same definition, $D(\omega_j) = \omega_0^2 - \omega_j^2 - 2i \omega_j \gamma$. Similarly, if we wanna get other response, we may substitute different frequencies, like for $\tilde{x}^{(2)}= x^{(2)}(2\omega_{1})e^{-2i \omega_{1}t}$, we can write $\frac{\mathrm{d}^2 \tilde{x}^{(2)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(2)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(2)} = \frac{-a \left( \frac{e E_2}{m} \right)^2 e^{-2i \omega_2 t}}{D^2(\omega_2)}$ The solution would be $x^{(2)}(2\omega_2) = \frac{-a \left( \frac{e}{m} \right)^2 E_2^2}{D(2\omega_2) D^2(\omega_2)}$ And for $\tilde{x}^{(2)}= x^{(2)}(\omega_{1}+\omega_{2})e^{-i (\omega_{1}+\omega_{2})t}$, $\frac{\mathrm{d}^2 \tilde{x}^{(2)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{x}^{(2)}}{\mathrm{d}t} + \omega_0^2 \tilde{x}^{(2)} = \frac{-a \left( \frac{e E_1}{m} \right)\left( \frac{e E_2}{m} \right) e^{-i (\omega_1+\omega_{2}) t}}{D(\omega_{1})D(\omega_2)}$ The solution is $x^{(2)}(\omega_1 + \omega_2) = \frac{-2a \left( \frac{e}{m} \right)^2 E_1 E_2}{D(\omega_1 + \omega_2) D(\omega_1) D(\omega_2)}$ All possible frequencies are $2\omega_{1}$, $2\omega_{2}$, $\omega_{1}+\omega_{2}$, $\omega_{1}-\omega_{2}$, $0$, and their opposite values. It is worth to note that negative frequencies correspond to complex conjugate of $E$. As a summary, the overall results are, $\begin{aligned} x^{(2)}(2\omega_2) &= \frac{-a \left( \frac{e}{m} \right)^2 E_2^2}{D(2\omega_2) D^2(\omega_2)}, \\ x^{(2)}(\omega_1 + \omega_2) &= \frac{-2a \left( \frac{e}{m} \right)^2 E_1 E_2}{D(\omega_1 + \omega_2) D(\omega_1) D(\omega_2)}, \\ x^{(2)}(\omega_1 - \omega_2) &= \frac{-2a \left( \frac{e}{m} \right)^2 E_1 E_2^*}{D(\omega_1 - \omega_2) D(\omega_1) D(-\omega_2)}, \\ x^{(2)}(0) &= \frac{-2a \left( \frac{e}{m} \right)^2 E_1 E_1^*}{D(0) D(\omega_1) D(-\omega_1)} + \frac{-2a \left( \frac{e}{m} \right)^2 E_2 E_2^*}{D(0) D(\omega_2) D(-\omega_2)}. \end{aligned}$ #### From displacement to susceptibility Recall the linear susceptibility, we have the definition $P^{(1)}(\omega_j) = \varepsilon_0 \chi^{(1)}(\omega_j) E(\omega_j)$ and for each electron, local polarization, $p^{(1)}(\omega_j) = -e x^{(1)}(\omega_j)$ So overall polarization, $P^{(1)}(\omega_j) = Np^{(1)}(\omega_j) = -N e x^{(1)}(\omega_j)$ and linear susceptibility, $\chi^{(1)}(\omega_j) = \frac{N e^2 / (\varepsilon_0 m)}{D(\omega_j)} = \frac{N e^2 / (\varepsilon_0 m)}{\omega_0^2 - \omega_j^2 - 2i \omega_j \gamma}$ If make the approximation for $D$, $D(\omega_j) = 2 \omega_0 (\omega_j - \omega_0 - i \gamma)$ we will get a simplified form and having real and imaginary parts separated. $\chi^{(1)}(\omega_j) = \frac{N e^2 / (2 \varepsilon_0 m \omega_0)}{\omega_j - \omega_0 - i \gamma} = \frac{N e^2}{2 \varepsilon_0 m \omega_0} \frac{(\omega_j - \omega_0) + i \gamma}{(\omega_j - \omega_0)^2 + \gamma^2}$ The real part is $\chi'(\omega_{j}) = \frac{N e^2}{2 \varepsilon_0 m \omega_0} \frac{(\omega_j - \omega_0)}{(\omega_j - \omega_0)^2 + \gamma^2}$ and imaginary part is $\chi''(\omega_{j})=\frac{N e^2}{2 \varepsilon_0 m \omega_0} \frac{i \gamma}{(\omega_j - \omega_0)^2 + \gamma^2}$. By writing the susceptibility in this way, the imaginary part, or absorption is exactly the Lorentzian form. (This is not the same as what we had in [[Classical Lorentz model#Polarization and relative permittivity under Lorentz model|previous version]], we did not make approximation) ![[Drawing 2024-09-09 21.52.43.excalidraw.svg]] The plot shows the normalization $\chi^{(1)}$. Real part and imaginary part are shown separately. >[!Note] >Recall the Lorentzian function follows >$\frac{A}{\pi} \left( \frac{\frac{\Gamma}{2}}{(x - x_0)^2 + \left( \frac{\Gamma}{2} \right)^2} \right)$ >Amplitude $A$ scales intensity; $x_{0}$ is the peak center; $\frac{\Gamma}{2}$ is HWHM (namely half of FWHM $\Gamma$). Here selecting $\frac{\Gamma}{2}=\gamma$ gives FWHM $2\gamma$. Do the same for second order susceptibility, for example, for $2\omega$, we now have the definition $P^{(2)}(2\omega_1) = \varepsilon_0 \chi^{(2)}(2\omega_1, \omega_1, \omega_1) E(\omega_1)^2$ and the calculated polarization $P^{(2)}(2\omega_1) = -N e x^{(2)}(2\omega_1)$ so the susceptibility, $\chi^{(2)}(2\omega_1, \omega_1, \omega_1) = \frac{N (e^3 / m^2) a}{\varepsilon_0 D(2\omega_1) D^2(\omega_1)}$ Notice that it can be written as the product of three linear susceptibility, $\chi^{(2)}(2\omega_1, \omega_1, \omega_1) = \frac{\varepsilon_0^2 m a}{N^2 e^3} \chi^{(1)}(2\omega_1) \left[ \chi^{(1)}(\omega_1) \right]^2$ This is very nice and we can estimation the magnitude of $\chi^{(2)}$, also, this nature is the reason why people write it as a function of three frequencies, like what we did in [[Nonlinear susceptibility]]. Check the definition in [[Understanding parametric nonlinear optical processes]]. For SFG, the definition of polarization is a little different with a numerical factor $2$, indicating distinct input field (for OR not the case, see the [[Understanding parametric nonlinear optical processes#^83ffdb|callout]]). This is the same for DFG and optical rectification. By combining $P^{(2)}(\omega_1 + \omega_2) = 2\varepsilon_0 \chi^{(2)}(\omega_1 + \omega_2, \omega_1, \omega_2) E(\omega_1) E(\omega_2)$ $P^{(2)}(\omega_1 + \omega_2) = -N e x^{(2)}(\omega_1 + \omega_2)$ we may get the susceptibility $\chi^{(2)}(\omega_1 + \omega_2, \omega_1, \omega_2))= \frac{N (e^3 / m^2) a}{\varepsilon_0 D(\omega_1 + \omega_2) D(\omega_1) D(\omega_2)} = \frac{\varepsilon_0^2 m a}{N^2 e^3} \chi^{(1)}(\omega_1 + \omega_2) \chi^{(1)}(\omega_1) \chi^{(1)}(\omega_2)$ Similarly, for DFG and OR, we have $\chi^{(2)}(\omega_1 - \omega_2, \omega_1, -\omega_2) = \frac{N (e^3 / m^2) a}{\varepsilon_0 D(\omega_1 - \omega_2) D(\omega_1) D(-\omega_2)} = \frac{\varepsilon_0^2 m a}{N^2 e^3} \chi^{(1)}(\omega_1 - \omega_2) \chi^{(1)}(\omega_1) \chi^{(1)}(-\omega_2)$ $\chi^{(2)}(0, \omega_1, -\omega_1) = \frac{N (e^3 / m^2) a}{\varepsilon_0 D(0) D(\omega_1) D(-\omega_1)} = \frac{\varepsilon_0^2 m a}{N^2 e^3} \chi^{(1)}(0) \chi^{(1)}(\omega_1) \chi^{(1)}(-\omega_1)$ With similar manner, one may continue the calculation and obtain 3rd and higher order susceptibility. >[!Note] >Miller's rule ([[OPTICAL SECOND HARMONIC GENERATION IN PIEZOELECTRIC CRYSTALS|from this paper]]) showed that for all non-centrosymmetric crystals >$\frac{\chi^{(2)}(\omega_1 + \omega_2, \omega_1, \omega_2)}{\chi^{(1)}(\omega_1 + \omega_2) \chi^{(1)}(\omega_1) \chi^{(1)}(\omega_2)}$ >are nearly constant. Compare with our previous result, we can see that if >$\frac{m a \varepsilon_0^2}{N^2 e^3}$ >is a constant, the given quantity would be constant. And for almost all atoms, the atomic number density $N$ is nearly the same, and other parameters except $a$ are constants. This nonlinear parameter $a$ can be estimated by the electron origin we used before. Here assumes linear and nonlinear displacement are in the same range, and the displacement could be estimated as the lattice constant $d$. Then we have >$m\omega_{0}^{2}d=mad^{2}$ >Then >$a=\frac{\omega_{0}^{2}}{d}$ >This gives an estimation of $\chi^{2}$ based on physical constants and lattice parameter, which is >$\chi^{(2)} = \frac{e^3}{\varepsilon_0 m^2 \omega_0^4 d^4}$ >And a numerical estimation would be, based on $d=0.3\ nm$ and $\omega_{0}=1\times 10^{16}\ rad/s$, $6.9\times 10^{-12}\ \frac{m}{V}$. ^bf9662 ### Centrosymmetric media #### System setup Now we may go one step further to see the origin of 3rd order nonlinear susceptibility. In theory, we still cannot make statement on its tensor property. But here we may assume that the material is isotropic, these type of centrosymmetric materials are quite common, like glasses. With this assumption, we may apply a vector field. And we assume we may neglect higher order potential functions besides what is defined, this also means that the electron displacement is relatively small. A simple form of restoring force for non-parabolic centrosymmetric potential would be $\tilde{\mathbf{F}}_{\text{restoring}} = -m\omega_0^2 \tilde{\mathbf{x}} + mb (\tilde{\mathbf{x}} \cdot \tilde{\mathbf{x}}) \tilde{\mathbf{x}}$ This $b$ is a parameter like $a$ we had before, a characterization for nonlinearity, and the potential energy (consider scalar form) would be $U(\tilde{x}) = - \int \tilde{{F}}_{\text{restoring}} d\tilde{x} = \frac{1}{2} m \omega_0^2 \tilde{x}^2 - \frac{1}{4} mb \tilde{x}^4$ ![[Drawing 2024-09-10 01.02.39.excalidraw.svg]] The plot shows how this potential looks like. >[!Notice] >In the restoring force expression, $mb (\tilde{\mathbf{x}} \cdot \tilde{\mathbf{x}}) \tilde{\mathbf{x}}$, the $\tilde{\mathbf{x}} \cdot \tilde{\mathbf{x}}$ gives a scalar and ensure the value, and last $\tilde{\mathbf{x}}$ ensures the direction. This is the same for later expression for optical fields. Write the Newton's second law, $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}} - b (\tilde{\mathbf{r}} \cdot \tilde{\mathbf{r}}) \tilde{\mathbf{r}} = -\frac{e \tilde{\mathbf{E}}(t)}{m}$ and set the field, $\tilde{\mathbf{E}}(t) = \mathbf{E}_1 e^{-i\omega_1 t} + \mathbf{E}_2 e^{-i\omega_2 t} + \mathbf{E}_3 e^{-i\omega_3 t} + \text{c.c.}$ or in the summation form $\tilde{\mathbf{E}}(t) = \sum_n \mathbf{E}(\omega_n) e^{-i\omega_n t}$. Here assume we have three distinct field cause we are considering 3rd nonlinearity. #### Solve the problem with perturbation ##### First order Just like before, write the perturbation parameter $\lambda$ and expand displacement $\tilde{\mathbf{r}}$. The problem becomes $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}} - b (\tilde{\mathbf{r}} \cdot \tilde{\mathbf{r}}) \tilde{\mathbf{r}} = -\lambda \frac{e \tilde{\mathbf{E}}(t)}{m}$ and $\tilde{\mathbf{r}}(t) = \lambda \tilde{\mathbf{r}}^{(1)}(t) + \lambda^2 \tilde{\mathbf{r}}^{(2)}(t) + \lambda^3 \tilde{\mathbf{r}}^{(3)}(t) + \cdots$ Insert this expanded $\tilde{\mathbf{r}}(t)$, and organize the equations by the order of $\lambda$, we should have the following equations. $\lambda$: $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}^{(1)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}^{(1)}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}}^{(1)} = -\frac{e \tilde{\mathbf{E}}(t)}{m}$ $\lambda^{2}$: $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}^{(2)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}^{(2)}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}}^{(2)} = 0$ $\lambda^{3}$: $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}^{(3)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}^{(3)}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}}^{(3)} - b (\tilde{\mathbf{r}}^{(1)} \cdot \tilde{\mathbf{r}}^{(1)}) \tilde{\mathbf{r}}^{(1)} = 0$ $\cdots$ Just like before, we can get the [[Classical Lorentz model]] result, just this time we have the vector field, and the displacement is the summation of all three frequencies, $\tilde{\mathbf{r}}^{(1)}(t) = \sum_n \mathbf{r}^{(1)}(\omega_n) e^{-i \omega_n t}$ where $\mathbf{r}^{(1)}(\omega_n) = \frac{-e \mathbf{E}(\omega_n) / m}{D(\omega_n)}$. Here the $D$ is still the complex denominator function, $D(\omega_n) = \omega_0^2 - \omega_n^2 - 2i \omega_n \gamma$ If we write them explicitly it, it will look like $\begin{aligned} \tilde{\mathbf{r}}^{(1)}(t) &= \sum_{n=1}^{3} \left( \frac{-e \tilde{\mathbf{E}}(\omega_n)}{m(\omega_0^2 - \omega_n^2 - 2i\gamma \omega_n)} e^{-i \omega_n t} \right) + \text{c.c.} \\ &= \left(\frac{-e \tilde{\mathbf{E}}(\omega_1)}{m(\omega_0^2 - \omega_1^2 - 2i\gamma \omega_1)} e^{-i \omega_1 t} + \frac{-e \tilde{\mathbf{E}}(\omega_2)}{m(\omega_0^2 - \omega_2^2 - 2i\gamma \omega_2)} e^{-i \omega_2 t} + \frac{-e \tilde{\mathbf{E}}(\omega_3)}{m(\omega_0^2 - \omega_3^2 - 2i\gamma \omega_3)} e^{-i \omega_3 t}\right) + \text{c.c.} \end{aligned} $ From the definition of polarization, we have $\mathbf{P}^{(1)}(\omega_n) = -N e \mathbf{r}^{(1)}(\omega_n)$ Treat different Cartesian components of polarization separately, $P_i^{(1)}(\omega_n)=-Ne r_{i}^{(1)}(\omega_{n})$ $P_i^{(1)}(\omega_n) = \varepsilon \sum_j \chi_{ij}^{(1)}(\omega_n) E_j(\omega_n)$ where $\chi_{ij}^{(1)}(\omega_n) = \chi^{(1)}(\omega_n) \delta_{ij}$, and $\chi^{(1)}(\omega_n) = \frac{N e^2 / m}{\varepsilon_0 D(\omega_n)}$. Here $\delta_{ij}$ is Kronecker delta. These $i,\ j$ are Cartesian components. If write explicitly, we will have $P_i^{(1)}(\omega_n) = \varepsilon_0 \left[ \chi_{ix}^{(1)}(\omega_n) E_x(\omega_n) + \chi_{iy}^{(1)}(\omega_n) E_y(\omega_n) + \chi_{iz}^{(1)}(\omega_n) E_z(\omega_n) \right],$ this is still for one input frequency $\omega_{n}$ (cause it is the first order, or one may consider it as a function of frequency, and it is). Here we see each component of the $\chi^{(1)}$ tensor, or $\chi^{(1)}(\omega_n) = \begin{pmatrix} \chi^{(1)}(\omega_n) & 0 & 0 \\ 0 & \chi^{(1)}(\omega_n) & 0 \\ 0 & 0 & \chi^{(1)}(\omega_n) \end{pmatrix}$ ##### Second order Check the second order expression, the right hand side is zero, indicating it's an equation with damping but no driven. The steady-state solution vanishes and $\tilde{\mathbf{r}}^{(2)} = 0$ ##### Third order Check the third order expression, for $\lambda^{3}$, the left and right hand side are $\frac{\mathrm{d}^2 \tilde{\mathbf{r}}^{(3)}}{\mathrm{d}t^2} + 2\gamma \frac{\mathrm{d}\tilde{\mathbf{r}}^{(3)}}{\mathrm{d}t} + \omega_0^2 \tilde{\mathbf{r}}^{(3)} = - \sum_{mnp} \frac{be^3 [\tilde{\mathbf{E}}(\omega_m) \cdot \tilde{\mathbf{E}}(\omega_n)] \tilde{\mathbf{E}}(\omega_p)}{m^3 D(\omega_m) D(\omega_n) D(\omega_p)} e^{-i(\omega_m + \omega_n + \omega_p)t}$ Define $\omega_{q}=\omega_m + \omega_n + \omega_p$ as the output frequency, the solution of above equation should have the form $\tilde{\mathbf{r}}^{(3)}(t) = \sum_q \mathbf{r}^{(3)}(\omega_q) e^{-i \omega_q t}$ Cause we only have three distinct input frequencies, so this $\sum_{q}$ does not have any significance. But things could be different if the distinct input fields are more than three. Recall the notation $(mnp)$ means the summation of the frequencies matches the output, i.e., $\omega_m + \omega_n + \omega_p$. So we just substitute $\tilde{\mathbf{r}}^{(3)}(t) = \sum_q \mathbf{r}^{(3)}(\omega_q) e^{-i \omega_q t}$ into the equation we wanna solve, we get $\left( -\omega_q^2 - i \omega_q 2 \gamma + \omega_0^2 \right) \mathbf{r}^{(3)}(\omega_q) = - \sum_{(mnp)} \frac{be^3 [\tilde{\mathbf{E}}(\omega_m) \cdot \tilde{\mathbf{E}}(\omega_n)] \tilde{\mathbf{E}}(\omega_p)}{m^3 D(\omega_q) D(\omega_m) D(\omega_n) D(\omega_p)}$ The coefficient $\left( -\omega_q^2 - i \omega_q 2 \gamma + \omega_0^2 \right)$ is just $D(\omega_{q})$, so $\mathbf{r}^{(3)}(\omega_q) = - \sum_{(mnp)} \frac{be^3 [\tilde{\mathbf{E}}(\omega_m) \cdot \tilde{\mathbf{E}}(\omega_n)] \tilde{\mathbf{E}}(\omega_p)}{m^3 D(\omega_q) D(\omega_m) D(\omega_n) D(\omega_p)}$ #### From displacement to susceptibility We now have the displacement, now consider polarization, $\mathbf{P}^{(3)}(\omega_q) = -N e \mathbf{r}^{(3)}(\omega_q)$ Also treat each Cartesian component separately, define the polarization on output direction $i$ as, $P_i^{(3)}(\omega_q) = \varepsilon_0 \sum_{ijkl} \sum_{(mnp)} \chi_{ijkl}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_j(\omega_m) E_k(\omega_n) E_l(\omega_p)$ > [!Note] > If expand this $P_{i}$, even if we do not consider the permutation of $mnp$ and the vector operation of the field, we will have 27 terms, namely > $P_i^{(3)}(\omega_q) = \varepsilon_0 \sum_{(mnp)} \begin{aligned} > & \chi_{ixxx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_x(\omega_n) E_x(\omega_p) > + \chi_{ixxy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_x(\omega_n) E_y(\omega_p) > + \chi_{ixxz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_x(\omega_n) E_z(\omega_p) \\ > + & \chi_{ixyx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_y(\omega_n) E_x(\omega_p) > + \chi_{ixyy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_y(\omega_n) E_y(\omega_p) > + \chi_{ixyz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_y(\omega_n) E_z(\omega_p) \\ > +& \chi_{ixzx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_z(\omega_n) E_x(\omega_p) > + \chi_{ixzy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_z(\omega_n) E_y(\omega_p) > + \chi_{ixzz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_x(\omega_m) E_z(\omega_n) E_z(\omega_p) \\ > +& \chi_{iyxx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_x(\omega_n) E_x(\omega_p) > + \chi_{iyxy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_x(\omega_n) E_y(\omega_p) > + \chi_{iyxz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_x(\omega_n) E_z(\omega_p) \\ > +& \chi_{iyyx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_y(\omega_n) E_x(\omega_p) > + \chi_{iyyy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_y(\omega_n) E_y(\omega_p) > + \chi_{iyyz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_y(\omega_n) E_z(\omega_p) \\ > +& \chi_{iyzx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_z(\omega_n) E_x(\omega_p) > + \chi_{iyzy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_z(\omega_n) E_y(\omega_p) > + \chi_{iyzz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_y(\omega_m) E_z(\omega_n) E_z(\omega_p) \\ > +& \chi_{izxx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_x(\omega_n) E_x(\omega_p) > + \chi_{izxy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_x(\omega_n) E_y(\omega_p) > + \chi_{izxz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_x(\omega_n) E_z(\omega_p) \\ > +& \chi_{izyx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_y(\omega_n) E_x(\omega_p) > + \chi_{izyy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_y(\omega_n) E_y(\omega_p) > + \chi_{izyz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_y(\omega_n) E_z(\omega_p) \\ > +& \chi_{izzx}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_z(\omega_n) E_x(\omega_p) > + \chi_{izzy}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_z(\omega_n) E_y(\omega_p) > + \chi_{izzz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) E_z(\omega_m) E_z(\omega_n) E_z(\omega_p) > \end{aligned} > $ > But since we have the dot product in the electric field, $\tilde{\mathbf{E}}(\omega_m) \cdot \tilde{\mathbf{E}}(\omega_n)$, perpendicular components will disappear. Therefore, we will have the first $\delta$, $\delta_{jk}$. > Then, because we need the output direction being $i$, therefore, the last field $\tilde{\mathbf{E}}(\omega_p)$ must point towards $i$. This gives the second $\delta_{il}$ > So one may see that a possible definition for susceptibility could be > $\chi_{ijkl}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p) = \frac{N b e^4 \delta_{jk} \delta_{il}}{\varepsilon_0 m^3 D(\omega_q) D(\omega_m) D(\omega_n) D(\omega_p)}$ > But this did not reflect the full symmetry due to the permutation nature. See later example. Notice than there are the summation over $mnp$, and the permutation of $m,\ n,\ p$ actually will repeat with the permutation of the Cartesian components. >[!Example] >For example, if we keep both summations, when permutation $m,\ n,\ p$ and $i,\ j,\ k,\ l$, we will encounter the following case: >$j=x,\ k=y,\ l=z$ >and the same time, we have >$\omega_{m}, \omega_{n}, \omega_{p}$ >So the fields will be $E_x(\omega_m) E_y(\omega_n) E_z(\omega_p)$, corresponds to $\chi_{ixyz}^{(3)}(\omega_q, \omega_m, \omega_n, \omega_p)$. >But after permutation, we will also have >$j=z,\ k=y,\ l=x$ >and the frequencies have the order >$\omega_{p}, \omega_{n}, \omega_{m}$ >This time we have the fields $E_z(\omega_p) E_y(\omega_n) E_x(\omega_m)$, this corresponds to $\chi_{izyx}^{(3)}(\omega_q, \omega_p, \omega_n, \omega_m)$. >Following our intuitions, changing the order of the same field should only be a semantic and notational difference, and will not change the physical quantity. And this is the intrinsic permutation. By accounting it, we will get the degeneracy coefficients in the susceptibility expression. It is conventional to define nonlinear susceptibilities in a manner that displays intrinsic permutation symmetry. Since there are six ($3!$) possible permutations of the order in which $E_j(\omega_m) E_l(\omega_n) E_k(\omega_p)$ may be taken, we define the $\frac{1}{6}$ of the sum of the six expressions of 3rd order susceptibility. It turns out that the final expression should be $\chi^{(3)}_{ijkl}(\omega_q, \omega_m, \omega_n, \omega_p) = \frac{N b e^4 \left[ \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right]}{3 \varepsilon_0 m^3 D(\omega_q) D(\omega_m) D(\omega_n) D(\omega_p)}$ This is a big tensor with 81 terms in total. However, it will be simplified in [[General symmetries and simplifications]]. We may also write it as a product of four linear susceptibility, namely $\chi^{(3)}_{ijkl}(\omega_q, \omega_m, \omega_n, \omega_p) = \frac{b m e^3 \varepsilon_0}{3 N^3 e^4} \left[ \chi^{(1)}(\omega_q) \chi^{(1)}(\omega_m) \chi^{(1)}(\omega_n) \chi^{(1)}(\omega_p) \right] \times \left[ \delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk} \right]$ And like what we did in [[#^bf9662|Miller's rule]], we can estimate the third order susceptibility by letting first and third order restoring force being at the same range, namely $m \omega_0^2 d = m b d^3$, we will have an estimation for $b = \frac{\omega_0^2}{d^2}$. If take $D(\omega)\approx\omega_{0}^{2}$, we have $\chi^{(3)} \simeq \frac{N b e^4}{\varepsilon_0 m^3 \omega_0^8} = \frac{e^4}{\varepsilon_0 m^3 \omega_0^6 d^5}$ Take $d=3\ nm$ and $\omega_{0}=7\times 10^{15}\ \frac{rad}{s}$, we have $\chi^{(3)} \simeq 344 \, \text{pm}^2/\text{V}^2$. This is an estimation given in [[Nonlinear Optics|page 32]].