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]:

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 $\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:

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.

 

Fig. 1 Dispersion relations of the present method, the red curve is the calculated result of Eq. (6) after normalization; the black curve is the numerical dispersion relations of free space.

Download Full Size | PDF

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₃ 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_i). 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:

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₁, k₂ are the wave vectors of the fundamental and SH wave respectively. The length and width of the OSL are 5mm and 3mm, and the focused point was designed to be the center of the exit surface. The expression $2k_{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 1064nm, 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.

 

Fig. 2 Numerical results with this differential method: (a) Domain structure of the LiTaO3 OSL; (b) Field distribution of the fundamental wave; (c) Field distribution of the second-harmonics.

Download Full Size | PDF

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].

 

Fig. 3 Intensity distribution of the harmonics at the focusing plane: (a) The fundamental wave; (b) Second-harmonic wave, respectively.

Download Full Size | PDF

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₃ 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.

 

Fig. 4 The field distributions of second-harmonic wave inside 23μm*23μm LiTaO₃ domain-inverted structure using FD method (a) and FDTD method (b), respectively. The wavelength was set to be 1064nm, and the working temperature was 100°C.The period of the structure was 7.7µm, the radius of cylinder was 1.93µm.

Download Full Size | PDF

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.