## Abstract

A finite difference method in real space is presented for calculating nonlinear optical processes in two-dimensional optical superlattices. The focused second-harmonic generation under the local quasi-phase-matched condition is calculated as an example. The field distribution of both the fundamental and the harmonic wave can be simulated well using this method, and the result agrees well with previous theoretical predictions and experimental studies. It is shown that this method is a simple and rapid technique to analysis nonlinear processes in optical superlattices.

©2012 Optical Society of America

## 1. Introduction

Since the theory of quasi-phase-matching (QPM) was introduced in 1962 [1], both theoretical and experimental studies on nonlinear processes in optical superlattices (OSL) have been investigated intensively [2–5]. In the early years, most of research in this field was focused on the one-dimensional (1D) case. In this situation, the nonlinear processes are basically governed by the 1D nonlinear coupled equations [6], which can be solved numerically without difficulties. Recently, studies on the nonlinear processes in two-dimensional (2D) [7–9] (or even three-dimensional (3D) [10]) optical superlattices have become a hot topic. The non-plane-wave interactions, such as the QPM processes involving Gaussian beams or Airy beams, have also attracted considerable attention [11–13]. Therefore, it is necessary and important to study the propagation and diffraction processes in a 2D optical superlattice. However, the 1D nonlinear coupled equations fail in this situation. It is well known that a lot of methods have been proposed to study the light transmission in photonic crystals for linear and nonlinear optical effects. The calculated techniques include the plane wave expansion method [14,15], the finite difference time domain (FDTD) method and the finite element method (FEM) [16–20]. FDTD and FEM can also be used to simulate second-order nonlinear effects in 2D systems. Among these methods, there are few reports on the calculation of the nonlinear processes in 2D or 3D optical superlattices, due to the long sample length (about several centimeters). In this paper, we report a simple and rapid calculation method which is suitable to calculate the nonlinear processes in 2D optical superlattices. This is also a finite difference (FD) method, but it solves the 2D paraxial wave equation directly in real space, which is several magnitudes faster than the FDTD method. For example, the focused second harmonic generation (SHG) process under the local QPM condition is calculated in detail, and the result agrees well with previous theoretical predictions and experimental studies.

## 2. Method

In order to elucidate the above idea, we consider the SHG processes in a 2D optical superlattice. The input fundamental wave is propagating along the *x* axis. Under the slowly varying envelope approximation, the evolution of the harmonics can be described by the paraxial wave equations [6]:

*K*is the coupling coefficient;

*f(x, y)*is the structure function of the 2D optical superlattice;

*k*,

_{1}, k_{2}*A*,

_{1}*A*are the wave vectors and amplitudes of the fundamental wave and second-harmonic wave respectively; is the phase mismatch between the fundamental and harmonic wave.

_{2}Along the *x* direction of the crystal, the propagating area is discretized by a grid mesh. Based on the Eqs. (1), the related difference equations can be written as:

Here $\mathrm{max}[\frac{\Delta x}{{k}_{1}\Delta {y}^{2}},\frac{\Delta x}{{k}_{2}\Delta {y}^{2}}]<\frac{1}{2}$, which means that if the condition is not satisfied, the iterative results will be divergent.

For given nonlinear materials, the coupling coefficient and the index of refraction at each grid point are known. Starting from appropriate initial conditions, the evolution of the fundamental wave and the second harmonics can be obtained by the above difference equations. When the input fundamental wave is a Gaussian beam, the corresponding initial conditions can be written as:

*A*,

_{1}*A*are the field values of the fundamental and the second harmonic respectively;

_{2}*A*is a constant;

*y*is the center of the fundamental wave and

_{0}*r*is the waist radius of the beam. The boundary condition for our method is also quite simple, which can be directly set to be zero:

_{0}*y = 0*and

*y = L*are the boundaries in

*y*-direction of the calculation area.

It should be note that Eq. (1) (the paraxial wave equation) is only valid for paraxial wave interactions, that is to say, the angle between the second harmonic (SH) wave and the *x*-axis should be small. This error is reflected in the dispersion relation of the differential equation, where the dispersion relation of Eq. (1) is:

Here ${k}_{ix}$and ${k}_{iy}$are the two components of along the x and y direction, and (i = 1,2) is the wave vector of the fundamental wave and SH wave respectively.

The discretization process in Eq. (2)-(3) will also bring numerical errors to the calculation; it can be also reflected in the dispersion relation. The numerical dispersion relation of Eq. (2)-(3) can be written as:

The dispersion relation of Eq. (1) is shown in Fig. 1
. The red curve is the calculated result of the Eq. (6), while the black curve is the numerical dispersion relations of free space. The curves corresponding to the fundamental and SH wave are exact same after normalization. It can be easily seen that the two curves almost coincide when the angle between the direction of the SH beam and the *x* axis is small, which means that the calculation is nearly accurate under this condition. Actually it was calculated that the angel should be smaller than 20 degrees.

## 3. Numerical results & discussion

In this section, we take the focused SHG process under the local QPM configuration as an example, realized in a 2D LiTaO_{3} OSL. Based on the Huygens-Fresnel principle for nonlinear interactions, the local QPM method is proposed by our group in 2008 [8], which is a general method for domain engineering with the conventional QPM method being a special case. It is considered that a small part of crystals can be regarded as a point source which emits the SH wave. The SH wave propagated from the point source (*x, y*) and SH beam can be focused at pre-designed point (*X _{i}*,

*Y*). To make the SH wave in phase at the focused point, the structure function (corresponding to the direction of domain polarization) of the optical superlattice to generate focused SHG can be determined:

_{i}Here $r=\sqrt{{({X}_{i}-x)}^{2}+{({Y}_{i}-y)}^{2}}$ is the distance between a point source and the focused SHG point, and *k _{1}*,

*k*are the wave vectors of the fundamental and SH wave respectively. The length and width of the OSL are 5

_{2}*mm*and 3

*mm*, and the focused point was designed to be the center of the exit surface. The expression $2{k}_{1}x+{k}_{2}r$ indicate the phase factor of SH wave, which determine the structure function and the direction of domain polarization inside the optical superlattice.

The domain structures of the optical superlattice are demonstrated in Fig. 2(a)
, which is obviously different from conventional 1D and 2D structures. Near the focused point, the domain boundary curves strongly, whereas the domain structure looks like the conventional far away from the focused point. The field distributions of the fundamental and the SH are shown in Fig. 2(b) and Fig. 2(c), respectively. The waist radius of the fundamental wave was about 300*μm*, and it propagated along the *x* axis of the OSL, with the incident point at the center of the OSL. The wavelength was set to be 1064*nm*, and the working temperature was 180°C. Figure 2(c) shows the field distribution of the SH wave, which was focused at the center of exit OSL surface. The propagation and diffraction properties are well resolved in the simulation. We can clearly see that the SH wave was focused gradually with the propagating of the fundamental wave and the focused point was at the exit surface of the OSL.

Moreover, the beam profiles of input fundamental wave and output SH wave at the focusing plane can also be calculated by current method. Figure 3 indicates the calculated results of the intensity of the fundamental and the SH at the focusing plane respectively. Be relative to the fundamental beam, the SH wave beam waist was considerable smaller. According to our calculation the waist of SH beam is about 40% of the fundamental beam and the conversion efficiency of SHG is about 21.4%, which agrees well with the theoretical predictions and experimental results in the previous work [8].

Now we would like to have a brief discussion on the advantages and disadvantages of the present method comparing to the conventional FDTD method. It is well known that the FDTD method is a powerful simulated technique, which is widely used to calculate the field distribution of propagating electromagnetic waves. Maxwell’s equations are subsequently discretized by using finite difference approximations in both spatial and temporal. While the present method used in this paper solves the paraxial wave equations directly, surely it is not appropriate for general 2D structures or structures with non-uniform refractive index. From the viewpoint of simplicity, it only needs to be discretized in real space, and the temporal complexity does not exist at all in this method.

On the other hand, in order to yield accurate results for a given wavelength, the grid spacing in FDTD method must be much less than the wavelength λ, typically λ/10 or smaller [21], while in present method the grid spacing can be set much larger (in our calculation and, which is even larger than the wavelength), which make it possible to simulate the field distribution of electromagnetic waves in the sample with several centimeters. Physically, under the varying envelope approximation, the counter propagating wave could be ignored, resulting in the allowable larger grid spacing. Also the calculation in general 2D structures can be simplified by using the paraxial wave equations. However, in FDTD method, the larger spacing step will lead to the divergence of calculation. Therefore, the current method is much simpler than the FDTD technique on the spatial complexity. Thus the method presented is a more efficient technique to calculate the 2D nonlinear process in OSL.

For detailed comparison, we simulate a 23μm*23μm LiTaO_{3} domain-inverted structure using present method and FDTD method, respectively. Actually it is almost impossible to calculate the structure shown in Fig. 2(a) using FDTD method, its computational time is enormous, reaching about 10 to the power of 8 seconds. Figure 4(a)-(b)
indicate the calculation results. The results are almost the same. However, the computational time in FD method is only one ten-thousandth of that in FDTD method.

This numerical calculation is a strict differential method based on the paraxial wave equations, thus the paraxial approximation is used to calculate the 2D nonlinear process. It cannot be applied to the general 2D structure. So a disadvantage of this method is that it can only be used when the angle of the fundamental and the SH wave is small (usually <20 degrees). Actually the error of the SH simulation will occur when the angle of the fundamental and the SH wave becomes bigger. At the angle of 20 degrees, the inaccuracy of wave vector is about two thousandth, which is still in the verge of tolerance.

## 4. Conclusion

In conclusion, a simple and rapid method was proposed to calculate the 2D nonlinear process using the centered difference method. Utilizing this method, the propagation and diffraction of the fundamental and harmonic waves can be simulated well. In principle, this technique is appropriate for calculating the nonlinear process and linear process in 2D or 3D OSLs. This method can be also applied to calculate other nonlinear processes such as sum-frequency generation and third-harmonic generation in optical superlattices. Besides, the method can also be extended to calculate the ultrafast pulse propagation in superlattices by combination of FD method and the split-step Fourier transform method.

## Acknowledgments

This work is supported by the State Key Program for Basic Research of China (No. 2010CB630703), the National Natural Science Foundation of China (No. 11074118 and No. 10874082) and the Doctoral Fund of Ministry of Education of China (No. 20090091110043).

## References and links

**1. **J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan, “Interactions between light waves in a nonlinear dielectric,” Phys. Rev. **127**(6), 1918–1939 (1962). [CrossRef]

**2. **A. Arie, G. Rosenman, V. Mahal, A. Skliar, M. Oron, M. Katz, and D. Eger, “Green and ultraviolet quasi-phase-matched second harmonic generation in bulk periodically-poled KTiOPO_{4},” Opt. Commun. **142**(4-6), 265–268 (1997). [CrossRef]

**3. **S. Wang, V. Pasiskevicius, F. Laurell, and H. Karlsson, “Ultraviolet generation by first-order frequency doubling in periodically poled KTiOPO_{4.},” Opt. Lett. **23**(24), 1883–1885 (1998). [CrossRef] [PubMed]

**4. **I. Yokohama, M. Asobe, A. Yokoo, H. Itoh, and T. Kaino, “All-optical switching by use of cascading of phase-matched sum-frequency generation and difference-frequency generation processes,” J. Opt. Soc. Am. B **14**(12), 3368 (1997). [CrossRef]

**5. **S. N. Zhu, Y. Y. Zhu, and N. B. Ming, “Quasi-phase-matched third-harmonic generation in a quasi-periodic optical superlattice,” Science **278**(5339), 843–846 (1997). [CrossRef]

**6. **R. W. Boyd, *Nonlinear Optics* (Elsevier Science, 2003).

**7. **V. Berger, “Nonlinear photonic crystals,” Phys. Rev. Lett. **81**(19), 4136–4139 (1998). [CrossRef]

**8. **Y. Q. Qin, C. Zhang, Y. Y. Zhu, X. P. Hu, and G. Zhao, “Wave-front engineering by Huygens-Fresnel principle for nonlinear optical interactions in domain engineered structures,” Phys. Rev. Lett. **100**(6), 063902 (2008). [CrossRef] [PubMed]

**9. **P. Xu, S. H. Ji, S. N. Zhu, X. Q. Yu, J. Sun, H. T. Wang, J. L. He, Y. Y. Zhu, and N. B. Ming, “Conical second harmonic generation in a two-dimensional χ^{(2)} photonic crystal: a hexagonally poled LiTaO_{3} crystal,” Phys. Rev. Lett. **93**(13), 133904 (2004). [CrossRef] [PubMed]

**10. **J. J. Chen and X. F. Chen, “Phase matching in three-dimensional nonlinear photonic crystals,” Phys. Rev. A **80**(1), 013801 (2009). [CrossRef]

**11. **C. Zhang, Y. Q. Qin, and Y. Y. Zhu, “Perfect quasi-phase matching for the third-harmonic generation using focused Gaussian beams,” Opt. Lett. **33**(7), 720–722 (2008). [CrossRef] [PubMed]

**12. **T. Ellenbogen, N. Voloch-Bloch, A. Ganany-Padowicz, and A. Arie, “Nonlinear generation and manipulation of Airy beams,” Nat. Photonics **3**, 395–398 (2009). [CrossRef]

**13. **I. Dolev, T. Ellenbogen, and A. Arie, “Switching the acceleration direction of Airy beams by a nonlinear optical process,” Opt. Lett. **35**(10), 1581–1583 (2010). [CrossRef] [PubMed]

**14. **M. S. Kushwaha, P. Halevi, G. Martínez, L. Dobrzynski, and B. Djafari-Rouhani, “Theory of acoustic band structure of periodic elastic composites,” Phys. Rev. B Condens. Matter **49**(4), 2313–2322 (1994). [CrossRef] [PubMed]

**15. **F. R. Montero de Espinosa, E. Jimenez, and M. Torres, “Ultrasonic band gap in a periodic two-dimensional composite,” Phys. Rev. Lett. **80**(6), 1208–1211 (1998). [CrossRef]

**16. **C. T. Chan, Q. L. Yu, and K. M. Ho, “Order-N spectral method for electromagnetic waves,” Phys. Rev. B Condens. Matter **51**(23), 16635–16642 (1995). [CrossRef] [PubMed]

**17. **Y. Tanaka, Y. Tomoyasu, and S. I. Tamura, “Band structure of acoustic waves in phononic lattices: Two-dimensional composites with large acoustic mismatch,” Phys. Rev. B **62**(11), 7387–7392 (2000). [CrossRef]

**18. **M. M. Sigalas and N. Garcia, “Importance of coupling between longitudinal and transverse components for the creation of acoustic band gaps: The aluminum in mercury case,” Appl. Phys. Lett. **76**(16), 2307 (2000). [CrossRef]

**19. **C. M. Reinke, A. Jafarpour, B. Momeni, M. Soltani, S. Khorasani, A. Adibi, Y. Xu, and R. K. Lee, “Nonlinear finite-difference time-domain method for the simulation of anisotropic,${\chi}^{(2)}$${\chi}^{(3)}$ ” J. Lightwave Technol. **24**(1), 624–634 (2006). [CrossRef]

**20. **A. Massaro, V. Tasco, M. T. Todaro, T. Stomeo, R. Cingolani, M. De Vittorio, and A. Passaseo, “FEM design and modeling of${\chi}^{(2)}$,” J. Lightwave Technol. **27**, 4262–4268 (2009). [CrossRef]

**21. **R. Drezek, A. Dunn, and R. Richards-Kortum, “A pulsed finite-difference time-domain (FDTD) method for calculating light scattering from biological cells over broad wavelength ranges,” Opt. Express **6**(7), 147–157 (2000). [CrossRef] [PubMed]