## Abstract

A versatile and accurate approach that combines a numerical iteration technique and a transfer-matrix method (TMM) is developed to solve the general problem of second harmonic generation (SHG) with pump depletion in quasi-phase-matched (QPM) nonlinear optical structures. We derive the iterative formulae from the nonlinear coupled wave equations and obtain the intensity distribution of fundamental wave and second harmonic wave by TMM. The approach shows quick numerical convergence of iteration and maintains perfect conservation of total energy. The simulation results show that the model coincides with the one under undepleted pump approximation very well when the SHG efficiency is small (well below 15%) and agrees very well with the effective nonlinear susceptibility model in handling general SHG problems even when the conversion efficiency is high up to 100%. Our method is applicable to general nonlinear optical structures, such as periodic, quasi-periodic, and aperiodic QPM structures, photonic crystals, and micro-cavities that might involve complicated modulation on the linear and nonlinear susceptibility.

© 2010 OSA

## 1. Introduction

In the context of nonlinear optical frequency conversion, particularly second-harmonic generation (SHG), the quasi-phase-matching (QPM) technique is known as an attractive way to obtain good phase matching and increase the nonlinear optical frequency conversion. This scheme was first proposed in 1962 [1] and has been experimentally investigated and confirmed in the one-dimensional (1D) or two-dimensional (2D), periodic or quasi-periodic structures [2–6]. Except for increasing the conversion efficiency, some other physical phenomena, such as compression of a tunable-chirp pulse [7], optical soliton [8], optical entanglement [9], the conical wave [10], enhanced Raman scattering [11] and optical focusing [12], have also been explored in the QPM structures. Recently, the generation and manipulation of Airy beams in an asymmetric quadratic nonlinear photonic crystal was reported [13]. To design proper QPM structures and analyze frequency conversion processes, a number of mathematic methods, such as the effective nonlinear susceptibility method under the undepleted pump approximation (ESMUP) [14,15], the Green function method [16] and the transfer matrix method (TMM) [17–21], have been developed. However, most of these methods are based on the undepleted pump approximation (UPA) and fail for the most interesting cases with high conversion efficiency. To circumvent the limitation, the effective nonlinear susceptibility method considering pump depletion (ESMPD) in the 1D periodic QPM structure was proposed [22]. According to the method, a homogeneous crystal with an effective nonlinear susceptibility can be used to model a QPM structure, where the ordinary analytical solution in the form of elliptical functions can yield very precise result compared with the exact one. The beam propagation method has been employed to calculate SHG in 2D and 3D nonlinear waveguide structures [23].

Other methods that start from the UPA solution and implement iterative algorithms to achieve gradual corrections to it were proposed to achieve better accuracy of solution in general SHG problems. In Ref [24,25]. a variational approach that combined finite-element method and fixed-point iteration was presented to discuss SHG in nonlinear optical films and 1D nonlinear photonic crystals, but the method still cannot handle the problem of very large SHG efficiency. More recently, a numerical iteration technique with the help of transfer-matrix method (TMM) was implemented to analyze 1D nonlinear photonic crystal micro-cavities [26]. The perturbation generated from pump depletion was obtained by assuming that the inhomogeneous term of the coupled-wave equation of fundamental wave (FW) was always small. However, when the conversion efficiency of second-harmonic wave (SHW) rises up to 50%, the inhomogeneous term becomes no longer small and the formulation become unstable and inaccurate in principle. Note that the SHG conversion efficiency reported by the paper is all below 15%, so the validity of the method in the high conversion efficiency cases still needs further confirmation.

In this paper, we propose a new iteration algorithm that is based on the TMM formulation that we developed in our previous literatures [19,20] to investigate SHG in the conventional 1D QPM structures. The algorithm has the advantage of fast numerical convergence, can maintain exact energy conservation, and more importantly, can successfully handle general structures with even high energy conversion. Numerical results show that our method coincides with the ESMUP when the conversion efficiency is below about 15% and agrees well with the ESMPD for general problems even when the conversion efficiency is very high up to 100% in QPM structures. In our model, pump depletion is taken into account and the usual slowly varying amplitude approximation (SVA) is not used. Therefore, it is accurate and applicable to a wide variety of one-dimensional nonlinear structures including the QPM ones, aperiodic structures, photonic crystals, and micro-cavities that might involve complicated modulation on the linear and nonlinear susceptibility. Moreover, the iteration algorithm can be extended to handle 2D and 3D nonlinear photonic crystals when it is incorporated with more general TMM formulation [21,27,28].

## 2. Theoretical model

Consider a usual 1D QPM structure consisting of a periodic array of alternating ferroelectric domain of positive and negative second-order nonlinear susceptibility. A plane wave FW is incident upon the structure. During the process of nonlinear optical interaction, a part of the FW energy is converted into SHW energy and leads to SHG in the forward and backward direction. Under the UPA, the solution of the nonlinear problem is straightforward by means of the conventional TMM [19,20]. Now consider a general situation where the pump depletion is no longer negligible and the conversion efficiency from FW to SHW might be high. The solution of the nonlinear problem is much more complicated. Benefiting from the idea of iteration and nonlinear TMM, we propose a kind of method to exactly solve the coupling and conversion problem of FW and SHW. The total sample of the 1D QPM structure can be divided into *N* slabs. The nonlinear coupled wave equations are given in the *n*th slab as follows [19,20,26],

In Eqs. (1) and (2) the physical quantities are defined as ${\beta}_{1n}={({k}_{1n}^{2}-{k}_{10}^{2}{\mathrm{sin}}^{2}\theta )}^{1/2}$, ${\beta}_{2n}={({k}_{2n}^{2}-{k}_{20}^{2}{\mathrm{sin}}^{2}\theta )}^{1/2}$, ${k}_{10}=\omega /c$, ${k}_{20}=2\omega /c$, ${k}_{1n}={n}_{1n}{k}_{10}$, and ${k}_{2n}={n}_{2n}{k}_{20}$. ${E}_{1n}$
$({E}_{2n})$and ${n}_{1n}$
$({n}_{2n})$ are respectively the electric field and refractive index for FW with angular frequency of *ω* (SHW with angular frequency of $2\omega $) in the *n*th slab, *c* is light speed in vacuum, *θ* represents the incident angle of FW, and${\chi}_{n}^{(2)}$ is the second-order coefficient (or nonlinear susceptibility) of the *n*th slab. The subscript * denotes the complex conjugate. The corresponding homogeneous equation of Eq. (1) is

We can easily obtain the general solution of Eq. (3). It includes the forward and backward terms of plane wave,

*n*th slice, and the sign + (−) denotes the forward-(backward-) propagating wave. By substituting Eq. (4) into Eq. (2) we can derive the solution of SHW as follows,

*n*th slab. Artificial infinitely-thin air films are brought in to surround each slab at the two ends in TMM [19,20], thus ${E}_{1n-1}^{\pm}$ $({E}_{2n-1}^{\pm})$ can be similarly defined as the plane wave coefficient of FW (SHW) in the air film at the left side of the

*n*th slab, and they are connected with ${\mathrm{\Omega}}_{1n}^{(\pm )}$(${\mathrm{\Omega}}_{2n}^{(\pm )}$) by using the boundary condition of electromagnetic field across the boundary. The transfer matrix for the

*n*th slice can be easily calculated by connecting ${E}_{1n-1}^{\pm}$ $({E}_{2n-1}^{\pm})$ with ${E}_{1n}^{\pm}$ $({E}_{2n}^{\pm})$. This single-slice transfer matrix is the basis to derive the transfer matrix for the unit cell and for the whole QPM structure. This is the basic principle of TMM for nonlinear optical structures [19,20]. As both the forward and backward propagating waves are considered in Eqs. (1) and (2), the TMM formulation is beyond the conventional slowly varying amplitude approximation [1].

Notice that Eqs. (4) and (5) serve as the basic formalisms of our nonlinear TMM under the undepleted pump approximation [19,20]. Here we call them the zero-order solution of iteration. In order to take into account the pump depletion, we substitute Eq. (4) and Eq. (5) back into Eq. (1), obtain the first-order solution of FW, and then derive the first-order solution of SHW from the first-order FW by substituting it into Eq. (2). This comprises a single run of iteration. In the iterative process, we only consider and reserve the $\mathrm{exp}[\pm i\mathrm{\Delta}{\beta}_{n}(z-{z}_{n-1})]$ terms and neglect the $\mathrm{exp}[\pm im\mathrm{\Delta}{\beta}_{n}(z-{z}_{n-1})]$
$(m=2,3,4\cdot \cdot \cdot )$ ones, assuming that they have little contribution because of strong phase mismatch involved in these terms. In fact, further consideration of these terms in our simulation almost does not make any difference of calculation results. After several iterations, the *i*th-order FW and SHW in the *n*th slab can be written as the following forms,

In Eqs. (6) and (7), $\mathrm{\Delta}{\beta}_{n}={\beta}_{2n}-2{\beta}_{1n}$ is the wave vector mismatch. The iteration starts from the zero-order solutions which are represented by the coefficients as,

The SHW coefficients in Eqs. (6) and (7) for the *i*th iteration are given by

In Eqs. (9) and (10) $i=0,1,2,3\cdot \cdot \cdot $, and the unknown parameters are defined as

*i*th-order solutions of FW and SHW in the

*n*th slab locating at ${z}_{n-1}$. Substituting the initial iteration conditions as expressed in Eq. (8) into Eqs. (6) and (7), we can find that the two expressions reduce to the zero-order solutions as shown in Eqs. (4) and (5). This means that the solutions expressed in Eqs. (6) and (7) are general and hold true for all iteration orders. Now, the iterative idea can be summarized as: first derive the

*i*th-order SHW from the

*i*th-order FW, and then obtain the (

*i +*1)th-order FW from the expressions of the

*i*th-order FW and SHW. Considering the initial iteration conditions, ${a}_{1n}^{(0)}={a}_{3n}^{(0)}={b}_{1n}^{(0)}={b}_{3n}^{(0)}=0$ and the initial boundary conditions, ${E}_{1N}^{(i)-}=0$, ${E}_{20}^{(i)+}=0$ and ${E}_{2N}^{(i)-}=0$, we describe the iterative process in detail as follows:

- (a) Find ${a}_{2n}^{(i)}$ and ${b}_{2n}^{(i)}$ $(i=0,1,2,3\cdot \cdot \cdot ,n=1,2,\cdot \cdot \cdot ,N)$ by TMM under the initial and continuous boundary conditions of FW;
- (b) Calculate ${g}_{1n}^{(i)}$, ${g}_{3n}^{(i)}$, ${f}_{1n}^{(i)}$, ${f}_{3n}^{(i)}$, ${h}_{1n}^{(i)}$, ${h}_{2n}^{(i)}$, and ${h}_{3n}^{(i)}$ $(i=0,1,2,3\cdot \cdot \cdot ,\text{}n=1,2,\cdot \cdot \cdot ,N)$ from Eqs. (9).1)-(9.7);
- (c) Find ${g}_{2n}^{(i)}$ and ${f}_{2n}^{(i)}$ $(i=0,1,2,3\cdot \cdot \cdot ,\text{}n=1,2,\cdot \cdot \cdot ,N)$ by TMM under the initial and continuous boundary conditions of SHW;
- (d) Calculate ${a}_{1n}^{(i+1)}$, ${a}_{3n}^{(i+1)}$, ${b}_{1n}^{(i+1)}$, and ${b}_{3n}^{(i+1)}$ $(i=0,1,2,3\cdot \cdot \cdot ,\text{}n=1,2,\cdot \cdot \cdot ,N)$ from Eqs. (10).1)- (10.4);
- (e) Replace
*i*by $i+1$ and repeat the process (a)-(d) until the convergence condition of Eq. (11) is satisfied.

*i +*1)th-order and

*i*th-order intensities of SHW at the

*z*-position of the

*n*th slab during the iteration.

*ε*is a sufficiently small number used to control the convergence of the iteration.

Both FW and SHW consist of forward and backward terms as expressed in Eqs. (6) and (7), so the forward and backward conversion efficiencies of SHW can be simultaneously calculated by

*f*and

*b*represents forward and backward, respectively.

During the process of frequency conversion, the energy conservation should always be obeyed. So when the convergence condition is satisfied in our model, we should have

If there is no reflectivity in the sample, for instance, when the QPM structure is connected with an index matched homogeneous medium, only the forward waves exist, and Eq. (14) becomes

## 3. Numerical results and discussion

To confirm the validity of our model, we investigate SHG in congruent lithium niobate (LiNbO_{3}, LN). The pump light of FW at the wavelength of 1.064 μm is normally incident upon the sample, and its field amplitude is chosen as 4.0 × 10^{6} V/m. The second-order coefficient of LN is set as 27.2pm/V [29], and the refractive index at room temperature can be derived as well [30].

Two optical structures are utilized to discuss our model: one is a homogeneous nonlinear LN film (Sample 1), and the other is a 1D periodic LN QPM structure (Sample 2). It is worth pointing out that, artificial infinitely-thin air films are brought in to surround each slab at the two ends in TMM [19,20], and the final results, such as conversion efficiency, are all expressed in the air background. Of course, if wanting to express the results in some other media such as LN, we can artificially introduce infinitely-thin films of the medium to surround each slab in TMM and assume that these films have no nonlinearity.

#### 3.1 Homogeneous optical film (Sample 1)

We first discuss the situation where the LN film is placed in a homogeneous medium background with a refractive index matched with that of LN where the reflectivity at the two ends of the sample becomes negligible. It is known that because LN has a very large dispersion between FW at 1.064 μm and SHW at 0.532 μm, where *n*
_{1} = 2.156 and *n*
_{2} = 2.2342, the sample is very much phase mismatched and has a small SHG conversion efficiency. In Fig. 1
, we plot the conversion efficiency as a function of the sample length or film thickness in order to compare our model with the UPA model and find that they become completely consistent with each other. This indicates that our model, like the UPA model, can also be implemented to analyze SHW in homogeneous optical films where the conversion efficiency is quite low.

We also investigate the LN film placed in the air background where the backward FW and SHW simultaneously take place, as is shown in Fig. 2 . Due to multi-reflection at the air and LN boundaries, many spikes appear in the curves of the forward and backward conversion efficiency of SHW versus the sample length, but the tendency of the forward SHW resembles the waveform in Fig. 1. Notice that the usual SVA model does not consider such a situation of impedance mismatch at the boundaries of nonlinear samples.

#### 3.2 Periodic QPM structure (Sample 2)

Sample 2 is an ideal periodic QPM structure with a duty circle of 50%, and each domain is precisely equal to the coherence length ($\pi /\mathrm{\Delta}\beta $~3.4016153 μm) so that the period is $2\pi /\mathrm{\Delta}\beta $ [14,15]. We have used our developed method to investigate SHG from such QPM structures. For comparison, the effective nonlinear susceptibility method under undepleted pump approximation (ESMUP) and the effective nonlinear susceptibility method considering pump depletion (ESMPD) are both employed to calculate the SHG conversion efficiency.

Figure 3
illustrates the conversion efficiency of SHW as a function of period number that is calculated by both the ESMUP and our model. Here the LN QPM structure is assumed to be index matched with the background medium. When the conversion efficiency of SHW is quite small or the sample length is quite short, our model coincides with the ESMUP very well. However, when the sample length gets longer, the two models begin to depart from each other when the conversion efficiency reaches a modest value of about 13.18%. When the length further increases, the conversion efficiency of SHW predicted in the ESMUP rises very quickly to go far beyond 100% at a large sample length. This is of course unphysical. As is well known, the UPA model is a successful one to reveal the basic physical picture of conversion efficiency in the QPM structure, such as the importance of reciprocal vector compensating phase mismatch in achieving high conversion efficiency. However, it becomes inapplicable to the high conversion efficiency case where the pump deletion should be considered. The curve of the conversion efficiency versus the sample length calculated by our model appears to follow a tanh^{2} shape which is well known in the case of phase-matched crystal [1].

We follow the authors of Ref [22]. and regard the periodic QPM structure as the phase-matched crystal whose effective nonlinear susceptibility is equal to the susceptibility of the same ideal crystal reduced by a factor of $2/\pi $. This is the ESMPD. Figure 4
displays the calculation result of the conversion efficiency versus period number by means of the ESMPD and our model. It is clear that our model keeps in good agreement with the ESMPD in a wide range of conversion frequency, no matter it is small or large. The curve predicted by our model is indeed a tanh^{2} function which is just the well-known analytical solution result. The conversion efficiency gets close to 100% for the long sample.

The above comparison study clearly indicates that our model is valid and reliable in handling general SHG problems in QPM structures. Note that the ESMPD is established in the periodic and quasi-periodic QPM structures where the nonlinear susceptibility can be written as the Fourier expansion series, so this simple model appears to be less efficient for aperiodic structures or random structures. In contrast, our model is not subject to such a difficulty and it can handle arbitrary structures with whatever complicated modulation of nonlinear susceptibility. This is because our model employs the TMM formulation, which in its nature deals with a structure slice by slice, slab by slab, whereas the slice or slab can have arbitrary thickness. For this reason, our model has a wider range of applicability than the UPA model, the ESMUP and ESMPD in handling general nonlinear optical structures, including periodic, quasi-periodic, and aperiodic ones.

We plot the efficiency of FW and SHW versus period number in Fig. 5
as well. The efficiency stands for transmission coefficient of FW and generation of SHW. From the figure, we can find that when the period number increases, the conversion efficiency of SHW gets growing while the transmission coefficient of FW falls down, but the total energy keeps invariable as an exact number of 1. The numerical calculation well confirms the prediction made by the model as in Eq. (15). Furthermore, when *N* = 2400, the conversion efficiency of SHW is 99.09% approaching unity, meaning that the total transfer of energy to SHW is possible in the current ideal QPM structure when the power of pump light is sufficiently large or when the sample is sufficiently long.

In Sample 2, the domain length is precisely set as the coherence length. Perfect phase matching is realized in the structure, so the conversion efficiency can reach a high value that is close to 100%. The situation will change if the domain length is set to have a little departure from the coherence length. The calculated conversion results for different domain lengths by means of our model are displayed in Fig. 6 . When the departure gets larger, the maximum conversion efficiency drops to a lower value. The curves resemble the well-known oscillation behavior that is predicted analytically for a phase-mismatched nonlinear crystal [1]. It is seen that our model can work well in this weakly phase-mismatched situation. The calculation results also support the general knowledge that the conversion efficiency of SHW is quite sensitive to the domain thickness in usual QPM structures.

We have investigated the conversion efficiency of SHG at different wavelengths of FW for a given QPM structure. The result is illustrated in Fig. 7 . We can see that the conversion efficiency of SHG vibrates with the incident wavelength, and at the wavelength of 1.064 μm, perfect phase matching is obtained and the conversion efficiency of SHW reaches the maximum. The bandwidth of SHG is quite small as about 0.1 nm. Calculation by means of the UPA for the small signal situation shows that the bandwidth has nearly the same small value. This small bandwidth is attributed to the strict request of phase matching between FW and SHW for a very long sample which is more than the incident wavelength by four orders of magnitude.

We consider another situation by placing Sample 2 in the air background instead of in the index-matched medium. The calculation results are plotted in Fig. 8 . The efficiencies of both the FW and SHW oscillate violently with respect to the sample length or period number, which is also induced by multiple-reflection at the two boundaries. As the current QPM structure is designed to compensate the phase mismatch between the forward FW and SHW, so the forward SHW plays a dominant role in the conversion process and its conversion efficiency is much larger than that of the backward SHW. Even with violent oscillation, the total energy of the four light waves still maintains perfect conservation, which again confirms Eq. (14) as predicted by our model. This is a further indication that our model works well in very complicated SHG problems.

## 4. Conclusion

In summary, we have developed a versatile and accurate approach that combines an efficient numerical iteration technique and the usual TMM formulation to solve the general problem of SHG in QPM nonlinear optical structures. Pump depletion of FW is taken into full account and the slowly varying amplitude approximation is not utilized in our model. The iteration algorithm exhibits quick numerical convergence and the conservation of total energy is maintained excellently even for very complicated QPM structures. Our numerical results show that our method coincides with the ESMUP very well when the conversion efficiency is below a value of about 15%, and agrees very well with the ESMPD in general situations even when the conversion efficiency is very large up to 100%. All of these features indicate that our model is effective, efficient, and accurate in handling SHG in general periodic QPM structures. As the method is based on the versatile and powerful TMM formulation, it is applicable to other general structures, such as aperiodic structures, photonic crystals, and micro-cavities that might involve complicated modulation on the linear and nonlinear susceptibility [31]. Moreover, it is expected that our model can be directly extended to handle nonlinear frequency conversion problems in 2D and 3D QPM and more complicated nonlinear optical structures.

## Acknowledgement

This work was supported by the National Natural Science Foundation at Nos. 10525419 and 10634080 and the State Key Development Program for Basic Research of China at No. 2007CB613205.

## 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. **G. D. Miller, R. G. Batchko, W. M. Tulloch, D. R. Weise, M. M. Fejer, and R. L. Byer, “42%-efficient single-pass cw second-harmonic generation in periodically poled lithium niobate,” Opt. Lett. **22**(24), 1834–1836 (1997). [CrossRef]

**3. **S. N. Zhu, Y. Y. Zhu, Y. Q. Qin, H. F. Wang, Ch. Zh. Ge, and N. B. Ming, “Experimental Realization of Second Harmonic Generation in a Fibonacci Optical Superlattice of LiTaO_{3},” Phys. Rev. Lett. **78**(14), 2752–2755 (1997). [CrossRef]

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

**5. **P. Ni, B. Ma, X. Wang, B. Cheng, and D. Zhang, “Second-harmonic generation in two-dimensional periodically poled lithium niobate using second-order quasiphase matching,” Appl. Phys. Lett. **82**(24), 4230–4232 (2003). [CrossRef]

**6. **B. Ma, T. Wang, Y. Sheng, P. Ni, Y. Wang, B. Cheng, and D. Zhang, “Quasiphase matched harmonic generation in a two-dimensional octagonal photonic superlattice,” Appl. Phys. Lett. **87**(25), 251103 (2005). [CrossRef]

**7. **A. M. Schober, G. Imeshev, and M. M. Fejer, “Tunable-chirp pulse compression in quasi-phase-matched second-harmonic generation,” Opt. Lett. **27**(13), 1129–1131 (2002). [CrossRef]

**8. **C. B. Clausen, O. Bang, and Y. S. Kivshar, “Spatial Solitons and Induced Kerr Effects in Quasi-Phase-Matched Quadratic Media,” Phys. Rev. Lett. **78**(25), 4749–4752 (1997). [CrossRef]

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

**10. **P. Xu, S. N. Zhu, X. Q. Yu, S. H. Ji, Z. D. Gao, G. Zhao, Y. Y. Zhu, and N. B. Ming, “Experimental studies of enhanced Raman scattering from a hexagonally poled LiTaO_{3} crystal,” Phys. Rev. B **72**(6), 064307 (2005). [CrossRef]

**11. **Y. Q. Qin, Ch. 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]

**12. **M. J. A. de Dood, W. T. M. Irvine, and D. Bouwmeester, “Nonlinear photonic crystals as a source of entangled photons,” Phys. Rev. Lett. **93**(4), 040504 (2004). [CrossRef]

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

**14. **J. D. McMullen, “Optical parametric interactions in isotropic materials using a phase-corrected stack of nonlinear dielectric plates,” J. Appl. Phys. **46**(7), 3076–3081 (1975). [CrossRef]

**15. **M. M. Fejer, G. A. Magel, D. H. Jundt, and R. L. Byer, “Quasi-phase-matched second harmonic generation: tuning and tolerances,” IEEE J. Quantum Electron. **28**(11), 2631–2654 (1992). [CrossRef]

**16. **X. H. Wang and B. Y. Gu, “Nonlinear frequency conversion in 2D χ^{(2)} photonic crystals and novel nonlinear double-circle construction,” Eur. Phys. J. B **24**(3), 323–326 (2001). [CrossRef]

**17. **D. S. Bethune, “Optical harmonic generation and mixing in multilayer media: analysis using optical transfer matrix techniques,” J. Opt. Soc. Am. B **6**(5), 910–916 (1989). [CrossRef]

**18. **Y. Jeong and B. Lee, “Matrix analysis for layered quasi-phase-matched media considering multiple reflection and pump wave depletion,” IEEE J. Quantum Electron. **35**(2), 162–178 (1999). [CrossRef]

**19. **J. J. Li, Z. Y. Li, and D. Z. Zhang, “Second harmonic generation in one-dimensional nonlinear photonic crystals solved by the transfer matrix method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **75**(5), 056606 (2007). [CrossRef]

**20. **J. J. Li, Z. Y. Li, Y. Sheng, and D. Z. Zhang, “Giant enhancement of second harmonic generation in poled ferroelectric crystals,” Appl. Phys. Lett. **91**(2), 022903 (2007). [CrossRef]

**21. **J. J. Li, Z. Y. Li, and D. Z. Zhang, “Nonlinear frequency conversion in two-dimensional nonlinear photonic crystals solved by a plane-wave-based transfer-matrix method,” Phys. Rev. B **77**(19), 195127 (2008). [CrossRef]

**22. **K. C. Rustagi, S. C. Mehendale, and S. Meenakshi, “Optical frequency conversion in quasi-phase-matched stacks of nonlinear crystals,” IEEE J. Quantum Electron. **18**(6), 1029–1041 (1982). [CrossRef]

**23. **H. M. Masoudi and J. M. Arnold, “Modeling second-order nonlinear effects in optical waveguides using a parallel-processing beam propagation method,” IEEE J. Quantum Electron. **31**(12), 2107–2113 (1995). [CrossRef]

**24. **G. Bao and D. C. Dobson, “Second harmonic generation in nonlinear optical films,” J. Math. Phys. **35**(4), 1622–1633 (1994). [CrossRef]

**25. **J. Yuan, “Computing for second harmonic generation in one-dimensional nonlinear photonic crystals,” Opt. Commun. **282**(13), 2628–2633 (2009). [CrossRef]

**26. **J. Xia, “Enhancement of second harmonic generation in one-dimensional nonlinear photonic-crystal microcavities,” Opt. Express **17**(22), 20069–20077 (2009). [CrossRef]

**27. **Z. Y. Li and L. L. Lin, “Photonic band structures solved by a plane-wave-based transfer-matrix method,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **67**(4), 046607 (2003). [CrossRef]

**28. **L. L. Lin, Z. Y. Li, and K. M. Ho, “Lattice symmetry applied in transfer-matrix methods for photonic crystals,” J. Appl. Phys. **94**(2), 811–821 (2003). [CrossRef]

**29. **V. G. Dmitriev, G. G. Gurazdyan, and D. N. Nikogosyan, *Handbook of nonlinear optical crystals* (Springer, Berlin, 1997), Vol. 64, p. 125.

**30. **G. J. Edwards and M. Lawrence, “A temperature-dependent dispersion equation for congruently grown lithium niobate,” Opt. Quantum Electron. **16**(4), 373–375 (1984). [CrossRef]

**31. **M. L. Ren and Z. Y. Li, “Giant enhancement of second harmonic generation in nonlinear photonic crystals with distributed Bragg reflector mirrors,” Opt. Express **17**(17), 14502–14510 (2009). [CrossRef]