Recently I have been working on the numerical calculation of second harmonic generation in nonlinear media. I refer to the book Laser Beam Propagation in Nonlinear Optical Media written by Shekhar Guha and Leonel P. Gonzalez, CRC Press 2013.
In chapter 7 of the book, they put forward a calculation algorithm to solve the equations below. $$\begin{align} \frac{\partial A_p}{\partial z} &= \frac{i}{2k_p} \nabla^2_{T}A_p + \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \frac{\partial A_s}{\partial z} &= \frac{i}{2k_s} \nabla^2_{T}A_s + \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$ Just simplify the above equations using linear and nonlinear operators $$\begin{align} \frac{\partial A_p}{\partial z} &= \hat{P_p}A_p + \hat{NL_p} \\ \frac{\partial A_s}{\partial z} &= \hat{P_s}A_s + \hat{NL_s} \end{align}$$ where the operators are $$\begin{align} \hat{P_p} &= \frac{i}{2k_p} \nabla^2_{T} \\ \hat{NL_p} &= \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \hat{P_s} &= \frac{i}{2k_s} \nabla^2_{T} \\ \hat{NL_s} &= \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$ Then it falls to the split step method to solve the equations.
Divide the propagation into N slices, each of length $\Delta z$
- Determine the incident pump field $A_p (x, y, z = 0, t)$
- Set $A_s(x,y,z=0,t) = 0$
- Set $\hat{NL} = 0$ and propagate the fields using $\hat{P}$ for a distance $\Delta z$. From $A_s(x,y,z=j\Delta z,t)$ to get $A_s(x,y,z=(j+1)\Delta z,t)$, we can use the the quasi-fast Hankel transform to get the results.
- Set $\hat{P} = 0$ and propagate the fields using $\hat{NL}$ for a distance $\Delta z$ using finite difference techniques (discussed below) to determine $A_P$ and $A_S$.
- The $A_P$ and $A_S$ fields become the inputs to the next $\Delta z$ slice.
- Repeat 3, 4 and 5 until fields propagate to the end of the crystal.
- Calculate the fields outside the crystal using the appropriate transmission coefficient.
Now comes the question:
Should I use the $A_s(x,y,(j+1)\Delta z,t)$ and $A_p(x,y,(j+1)\Delta z,t)$ that I calculated in step 3, or in step 4?
The detail of step 4 described in the book is listed here: Setting $\hat{P} = 0$, solve the nonlinear portion of the equations.
$$\begin{align} \frac{\partial A_p}{\partial z} &= \frac{i2d_\text{eff} \omega_p}{cn_p}A^*_p A_s e^{i(K_s - 2K_p)z} -\frac{\alpha_p}{2}A_p \\ \frac{\partial A_s}{\partial z} &= \frac{i2d_\text{eff} \omega_p}{cn_s}A_p A_p e^{-i(K_s - 2K_p)z} -\frac{\alpha_s}{2}A_p \end{align}$$
These equations are solved using the implicit and forward finite difference method from the known values expressed in a finite difference form as $$\begin{align} \frac{A_s^{(j+1)\Delta z} - A_s^{(j)\Delta z}}{\Delta z} &= C_s A^{(j+1)\Delta z}_p A_p^{(j+1)\Delta z} -\frac{\alpha_s}{2}A^{(j+1)\Delta z}_s \\ \frac{A_p^{(j+1)\Delta z} - A_p^{(j)\Delta z}}{\Delta z} &= C_s A^{*,j\Delta z}_p A_s^{(j+1)\Delta z} -\frac{\alpha_p}{2}A^{(j+1)\Delta z}_p \end{align}$$
where
$$\begin{align} C_p &= \frac{i2d_\text{eff} \omega_p}{cn_p} e^{i(K_s - 2K_p)(j \Delta z)} \\ C_s &= \frac{i2d_\text{eff} \omega_p}{cn_s} e^{-i(K_s - 2K_p)(j \Delta z)} \end{align}$$
It confused me that I have already got the calculated result of the $A_s(x,y,(j+1)\Delta z,t)$ and $A_p(x,y,(j+1)\Delta z,t)$ through step 3. And I have to calculate the $A_s(x,y,(j+1)\Delta z,t)$ and $A_p(x,y,(j+1)\Delta z,t)$ again in step 4 using the result of the $A_s(x,y,(j+1)\Delta z,t)$ and $A_p(x,y,(j+1)\Delta z,t)$ that already calculated in step 3.
Can anyone help me on this problem? What's the relationship between step 3 and step4?