## Sum frequency generation under wave equation description ### Simple case: SFG in lossless medium with collimated, monochromatic and continues-wave inputs Consider the most simple case. The waves are collimated, monochromatic and continues. The medium is lossless. And we have normal incidence only. This is shown in the following sketch. ![[Drawing 2024-09-24 00.57.51.excalidraw.svg]] Still adapt the time domain description, which means we should consider the equation $\nabla^2 \tilde{\mathbf{E}}_n - \frac{\boldsymbol{\varepsilon}^{(1)}(\omega_n)}{c^2} \cdot \frac{\partial^2 \tilde{\mathbf{E}}_n}{\partial t^2} = \frac{1}{\varepsilon_0 c^2} \frac{\partial^2 \tilde{\mathbf{P}}_{\mathrm{NL},n}}{\partial t^2}$ In this simple case, we can ignore the tensor nature and write all quantities as scalars. $\nabla^2 \tilde{{E}}_n - \frac{{\varepsilon}^{(1)}(\omega_n)}{c^2} \cdot \frac{\partial^2 \tilde{{E}}_n}{\partial t^2} = \frac{1}{\varepsilon_0 c^2} \frac{\partial^2 \tilde{{P}}_{\mathrm{NL},n}}{\partial t^2}$ As stated before in [[Nonlinear wave equation]], the nonlinear polarization acts as a source term. If we do not have this part, namely we are in linear medium, the solution for $\omega_{3}$ would just be like a free space field with the expression, $\tilde{E}_3(z,t) = A_3 e^{i(k_3 z - \omega_3 t)} + \text{c.c.}$ And the amplitude $A_{3}$ is a **constant**, the wavevector $k_{3}=\frac{n_{3}\omega_{3}}{c}$, refractive index $n_3^2 = \varepsilon^{(1)}(\omega_3)$. But now we have this nonlinear polarization, we should modified the expression of the field. When the nonlinear term is not so large, we expect to have the same sinusoidal form, but the amplitude $A_{3}$ becomes **a slow varying function of space $z$**, here we may consider it as a slow varying envelope. >[!Notice] >This amplitude or say envelope function has its own phase (this is intuitive since this is also a periodic function), which could differ from the the phase of the fast varying part (namely in $e^{i(k_3 z - \omega_3 t)}$). >This phase will be shown later during the derivation. Now we try to figure out the expression of $A_{3}$. #### Write the input field and nonlinear polarization Write the nonlinear polarization in above wave equation as $\tilde{P}_{3}$, and $\tilde{P}_3(z,t) = P_3 e^{-i\omega_3 t} + \text{c.c.}$ And from [[Understanding parametric nonlinear optical processes#Sum and difference frequency generation]], we see that $P(\omega_1 + \omega_2) = 2 \varepsilon_0 \chi^{(2)} E_1 E_2$, with the $d$ notation, it is $P_3 = 4 \varepsilon_0 d_{\text{eff}} E_1 E_2$ These $E_{1}$ and $E_{2}$ are spatially varying, namely, $\tilde{E}_i(z,t) = E_i e^{-i \omega_i t} + \text{c.c.}$ where $E_{i}=A_{i}e^{ik_{i}z}$. This $A_{i}$ is the amplitude of input field (or say, envelope), $i=1,2$. So $P_{3}$ is also spatially varying. This makes $P_{3}$ being $P_3 = 4 \varepsilon_0 d_{\text{eff}} A_1 A_2 e^{i(k_1 + k_2)z} \equiv p_3 e^{i(k_1 + k_2)z}$ Here this $p_{3}$ depends on the complex amplitudes of two input fields. The $\tilde{P}_3(z,t)$ then becomes $\tilde{P}_3(z,t)=4 \varepsilon_0 d_{\text{eff}} A_1 A_2 e^{i(k_1 + k_2)z}e^{-i\omega_3 t} + \text{c.c.} $ >[!Notice] >In the derivation we do not require $A_{1}$ and $A_{2}$ being $z$ independent. Although we removed the sinusoidal part, but this amplitude could still be varying, just like $A_{3}$. >We also do not guarantee that these amplitudes are real. This makes sense cause they are periodic and naturally has their own (spatial) phase, so using complex function would be intuitive. They are not simple plane waves and remember to consider these amplitude $A_{1}$, $A_{2}$ and $A_{3}$ as **envelope functions**. >For simplification, in many cases we may consider these $A_{1}$ and $A_{2}$ being constant. This is valid when the conversion of the input fields into the sum-frequency field is not too large. #### Insert into wave equation Now insert the expression $\tilde{P}_3(z,t)=4 \varepsilon_0 d_{\text{eff}} A_1 A_2 e^{i(k_1 + k_2)z}e^{-i\omega_3 t} + \text{c.c.}$ and electric field $\tilde{E}_3(z,t) = A_3 e^{i(k_3 z - \omega_3 t)} + \text{c.c.}$ into the wave equation, we have $\left[ \colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 A_3}{\mathrm{d}z^2}\)} + \colorbox{#ffa94d66}{\(2i k_3 \frac{\mathrm{d} A_3}{\mathrm{d}z}\)} - \colorbox{#4dff5062}{\(k_3^2 A_3 \)} + \colorbox{#ff89a962}{\(\frac{\varepsilon^{(1)}(\omega_3)\omega_3^2 A_3}{c^2}\)} \right] e^{i(k_3 z - \omega_3 t)} + \text{c.c.}= \frac{-4 d_{\text{eff}} \omega_3^2}{c^2} A_1 A_2 e^{i[(k_1 + k_2)z - \omega_3 t]} + \text{c.c.}$ The $\colorbox{#4dff5062}{third}$ and $\colorbox{#ff89a962}{forth}$ terms are cancelled. This is because $k_{3}=\frac{n_{3}\omega_{3}}{c}$, $k_{3}^{2}=\frac{n_{3}^{2}\omega_{3}^{2}}{c^{2}}=\frac{\varepsilon^{(1)}(\omega_3)\omega_{3}^{2}}{c^{2}}$. The equality still holds even if we remove the complex conjugate. >[!Notice] >For two random functions $f$ and $g$, and an equality $fe^{ik_{1}x}+\text{c.c.}=ge^{ik_{2}x}+\text{c.c.}$, removing the complex conjugate the equality will **NOT** automatically hold, even if $f$ and $g$ are real. The proof is trivial and will not be shown here. >However, here we may remove the complex conjugate is purely based on our physical assumptions. We are considering a specific SFG process, the output frequency is $\omega_{3}$ only and the propagating direction is only $k_{3}$. The input fields also have predetermined directions. Therefore, we could remove the negative frequency and negative wavevectors without violating the equality. $\left[ \colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 A_3}{\mathrm{d}z^2}\)} + \colorbox{#ffa94d66}{\(2i k_3 \frac{\mathrm{d} A_3}{\mathrm{d}z}\)}\right] e^{i(k_3 z - \omega_3 t)} = \frac{-4 d_{\text{eff}} \omega_3^2}{c^2} A_1 A_2 e^{i[(k_1 + k_2)z - \omega_3 t]}$ Remove $e^{-i\omega_{3}t}$, $ \colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 A_3}{\mathrm{d}z^2}\)} + \colorbox{#ffa94d66}{\(2i k_3 \frac{\mathrm{d} A_3}{\mathrm{d}z}\)} = \frac{-4 d_{\text{eff}} \omega_3^2}{c^2} A_1 A_2 e^{i(k_1 + k_2-k_3)z}$ Under [[Slowly varying envelope approximation]], we may ignore the $\colorbox{#a5d8ff66}{first}$ term because $\left| \colorbox{#a5d8ff66}{\(\frac{\mathrm{d}^2 A_3}{\mathrm{d}z^2}\)} \right| \ll \left| \colorbox{#ffa94d66}{\(2i k_3 \frac{\mathrm{d} A_3}{\mathrm{d}z}\)} \right|$ This holds when the variation of $A_{3}$ should be significantly smaller than $A_{3}$ itself in a wavelength $\lambda_{3}$. With this slow varying approximation applied, we have $\frac{\mathrm{d} A_3}{\mathrm{d}z} = \frac{2i d_{\text{eff}} \omega_3}{n_3 c} A_1 A_2 e^{i \Delta kz}$ $k_{3}$ is replaced with $k_{3}=\frac{n_{3}\omega_{3}}{c}$, and $\Delta k = k_1 + k_2 - k_3$, the wavevector mismatch. This equation is called coupled-amplitude equation, for it couples the amplitude of $\omega_{3}$ to that of $\omega_{1}$ and $\omega_{2}$. Note that since in this SFG process, all waves have to satisfy the wave equation. So we may follow the same derivation and calculate the amplitude of input fields, $A_{1}$ and $A_{2}$. Just this time we could have negative directions, and they looks like $\frac{\mathrm{d}A_1}{\mathrm{d}z} = \frac{2i d_{\text{eff}} \omega_1}{n_1 c} A_3 A_2^* e^{-i \Delta k z}$ $\frac{\mathrm{d}A_2}{\mathrm{d}z} = \frac{2i d_{\text{eff}} \omega_2}{n_2 c} A_3 A_1^* e^{-i \Delta k z}$ One simplification writing these equations are replace these complex amplitudes (envelopes) $A_{i}$ with the amplitude of polarization, $p_{i}$. For $p_{3}$, we had $P_3 = 4 \varepsilon_0 d_{\text{eff}} A_1 A_2 e^{i(k_1 + k_2)z} \equiv p_3 e^{i(k_1 + k_2)z}$, so $p_{3}=4 \varepsilon_0 d_{\text{eff}} A_1 A_2$. Similarly, $p_{1}=4 \varepsilon_0 d_{\text{eff}} A_3 A_2^{*}$, $p_{2}=4 \varepsilon_0 d_{\text{eff}} A_3 A_1^{*}$. These are also complex amplitudes (envelopes). And we may get something like $\frac{\mathrm{d}A_j}{\mathrm{d}z} = \frac{i \omega_j}{2 \varepsilon_0 n_j c} p_j e^{i \Delta k z}$ $j$ denotes for field $j$ with frequency $\omega_{j}$. Also we are considering lossless medium, so we may apply [[General symmetries and simplifications#Full permutation symmetry time reversal symmetry for lossless medium|full permutation symmetry]] and the coupling coefficient $d_{\text{eff}}$ has the same value in each equation. That's why we may use the same $d_{\text{eff}}$ in $p_{j}$ expressions. One may see here even if we have almost perfect setup (constant higher-order susceptibilities and extremely well defined boundary conditions), we still have coupled wave equations, and such complexity arises from the coupled PDEs themselves. #### Phase-matching considerations, amplitude and intensity Now we examine the wavevector mismatch $\Delta k = k_1 + k_2 - k_3$. For simplicity, assume input fields are not too large so they can be considered as **constants**. **Case 1**: perfect phase matching If we have perfect wavevector matching then $\Delta k = k_1 + k_2 - k_3=0$. From $\frac{\mathrm{d} A_3}{\mathrm{d}z} = \frac{2i d_{\text{eff}} \omega_3}{n_3 c} A_1 A_2 e^{i \Delta kz}$ we see that $\frac{\mathrm{d} A_3}{\mathrm{d}z} = \text{constant}$, and $A_{3}$ increases linear with $z$. Under this condition, all atomic dipoles are properly phased and perfect constructive interference is achieved. The generated wave maintain a fixed phase relation with the nonlinear polarization and the efficiency is maximized. **Case 2**: phase mismatch If perfect wavevector matching does not satisfy, we have to integrate along the material to get $A_{3}$. $\mathrm{d} A_3 = \frac{2i d_{\text{eff}} \omega_3}{n_3 c} A_1 A_2 e^{i \Delta kz} \mathrm{d}z$ $A_3(L) = \frac{2i d_{\text{eff}} \omega_3 A_1 A_2}{n_3 c} \int_0^L e^{i \Delta k z} \, \mathrm{d}z = \frac{2i d_{\text{eff}} \omega_3 A_1 A_2}{n_3 c} \left( \frac{e^{i \Delta k L} - 1}{i \Delta k} \right)$ Calculate the intensity, we have $I_3 = 2 n_3 \varepsilon_0 c |A_3|^2$ And for $I_{1}$, $I_{2}$, $I_1 = 2 n_1 \varepsilon_0 c |A_1|^2$ $I_2 = 2 n_2 \varepsilon_0 c |A_2|^2$ Insert the expression of $A_{3}$, and replace $|A_1|^2$, $|A_2|^2$ with corresponding intensities, $\begin{aligned} I_3 &= \frac{8 n_3 \varepsilon_0 d_{\text{eff}}^2 \omega_3^2 |A_1|^2 |A_2|^2}{n_3^2 c} \left| \frac{e^{i \Delta k L} - 1}{\Delta k} \right|^2\\ &= \frac{8 \varepsilon_0 d_{\text{eff}}^2 \omega_3^2 |A_1|^2 |A_2|^2}{n_3 c}L^2 \left( \frac{e^{i \Delta k L} - 1}{\Delta k L} \right)\left( \frac{e^{-i \Delta k L} - 1}{\Delta k L} \right)\\ & =\left(\frac{8\varepsilon_{0}d_{\text{eff}}^{2}\omega_{3}^{2}}{n_{3}c} \right) \left( \frac{I_{1}}{2n_{1}\varepsilon_{0}c}\right)\left( \frac{I_{2}}{2n_{2}\varepsilon_{0}c}\right) \left[ 2 L^2 \frac{1 - \cos(\Delta k L)}{(\Delta k L)^2}\right] \\ &= \frac{2 d_{\text{eff}}^2 \omega_3^2 I_1 I_2}{n_1 n_2 n_3 \varepsilon_0 c^3}L^2 \frac{\sin^2(\Delta k L / 2)}{(\Delta k / 2)^2} \\ &= \frac{2 d_{\text{eff}}^2 \omega_3^2 I_1 I_2}{n_1 n_2 n_3 \varepsilon_0 c^3} L^2 \, \text{sinc}^2 \left( \frac{\Delta k L}{2} \right) \end{aligned}$ >[!Note]- >Under our definition, the intensity is given by >$I_i = 2 n_i \varepsilon_0 c |A_i|^2, \quad i = 1, 2, 3$ >The intensity is given by the magnitude of the time-averaged Poynting vector. Sometimes people write >$\langle I \rangle = \frac{n c \epsilon_0}{2} \langle |E_0(t)|^2 \rangle$ >Where $E(t)=E_{0}(t)f(\omega t+\phi)$, $f$ is a real function. The difference is just because we use the complex form, in expressions above $\tilde{E}_i(z,t) = E_i e^{-i \omega_i t} + \text{c.c.}$. So the amplitude $A_{i}$ is half of $E_{0}$, if we write it as real function form. From the expression of $I_{3}$, the intensity of output wave, one can see that it only varies with the $\text{sinc}^2 \left( \frac{\Delta k L}{2} \right)$ part. The plot is shown below. ![[Drawing 2024-09-24 15.14.41.excalidraw.svg]] Only at wavevector matching condition we could obtain the highest intensity. With $|\Delta k|L$ increases, the efficiency decreases with some oscillations. A characteristic length $L_{\text{coh}}$, coherent buildup length, is defined, separating phase-matching and mismatching, while the latter means that the output wave can get out of phase with the driven polarization, which typically happens $L$ is larger than approximately $\frac{1}{\Delta k}$, and $L_{\text{coh}}=\frac{2}{\Delta k}$ Note this definition is purely out of mathematical convenience. With $L_{\text{coh}}$, the $\text{sinc}^2 \left( \frac{\Delta k L}{2} \right)$ becomes $\text{sinc}^2 \left(L / L_{\text{coh}} \right)$. The dashed line in above plot indicates the position of $\frac{L}{L_{\text{coh}}}=1$.