Optimal algorithms for the average-constrained maximum-sum segment problem

Quasi-3-D Beam-Propagation Method for Modeling

Nonlinear Wavelength Conversion

Shing Mou, Ching-Fuh Lin, Senior Member, IEEE, Member, OSA, and Hsu-Feng Chou

Abstract—The two-dimensional (2-D) iterative finite dif-ference beam-propagation method (IFD-BPM) is modified to model the cylindrically symmetric three-dimensional (quasi-3-D) second-order nonlinear wavelength conversion in quasi-phase-matched condition. The study shows that the dif-ference between the 2-D and 3-D schemes is small for the guided waves but large for the nonguided beams. The comparison with experimental results shows that the quasi-3-D IFD-BPM is closer to reality than the 2-D scheme. In addition, simulation using the quasi-3-D IFD-BPM reveals that plane-wave and Gaussian-beam assumptions are not sufficient for estimating the nonlinear con-version and beam propagation in second-order nonlinear devices. Index Terms—Difference frequency generation, finite-difference beam-propagation method, periodic poled lithium niobate, quasi-phase-matched techniques.

I. INTRODUCTION

Q

UASI-PHASE-MATCHED (QPM) frequency conversion [1], [2], using periodically poled ferroelectric crystals, has various applications because of its versatile and efficient conversion abilities. QPM devices based on LiNbO , LiTaO , and KTP [3]–[5] hold promise for the generation of a wide range of optical frequencies that are otherwise difficult to obtain. For example, the QPM second harmonic generation (SHG) [6] provides coherent light sources in the visible region [7]. Other QPM second-order nonlinear effects like sum frequency gen-eration (SFG) and difference frequency gengen-eration (DFG) are useful in frequency tripling and quadrupling [8], mid-infrared (IR) generation [9], [10], and wavelength-division multiplexing (WDM) [11]. Large nonlinear phase shifts induced by QPM second-order nonlinear effect [12] have also been proposed to realize all-optical switching [13] and short-pulse compression [14].

The advantage of QPM techniques over birefrin-gent-phase-matched techniques is the possibility to phase-match the full transparency range of the material and the availability of large nonlinear coefficients [1] through the engineerable periodic structure. To design and engineer those devices, it is important to understand the functions of these devices in advance, so theoretical models are important tools. However, analytical models usually require many approx-imations like the nondepletion of the pump wave or plane wave assumption [6], [15]–[18]. In addition, when large or irregular

Manuscript received March 3, 2000; revised December 29, 2000.

The authors are with the Department of Electrical Engineering and Graduate Institute of Electro-Optical Engineering, National Taiwan University, Taipei, Taiwan, R.O.C.

Publisher Item Identifier S 0733-8724(01)03638-6.

geometrical variations exist or the depletion of the pump wave is not negligible, precise analytical modeling of these devices is quite difficult. As a result, numerical methods become neces-sary for general and accurate analyses. The beam-propagation method (BPM) [19] is a powerful and flexible approach to design and simulate optical devices. Applications of BPM to linear devices had been studied extensively and extended to devices with second-order nonlinear effects, mostly SHG. The simulation is based on schemes like fast Fourier transform (FFT) [20] and finite-element (FE) [21]–[23], and finite-differ-ence (FD) methods [24]–[30]. FFT-BPM is preferred less for its low efficiency and accuracy. Its improvements are still under investigation [19], [31]. On the other hand, the comparison of FD-BPM and FE-BPM makes the former more attractive due to its simplicity of implementation.

To further improve the accuracy of FD-BPM in the case of nonlinear wavelength conversion, iterative finite-difference BPM (IFD-BPM) was proposed [29], [30]. IFD-BPM was shown to have good reliability, efficiency, and accuracy in com-parison with other FD-BPMs. Such IFD-BPM was developed in two-dimensional (2-D) cases before for lucid manifestation of comparisons among those BPM algorithms. However, the 2-D model does not take into account the beam divergence in both transverse directions. When the propagation beam is not confined in the waveguide, the 2-D model overestimates the nonlinear conversion. In order to make the simulation close to the real situation in general, three-dimensional (3-D) consid-eration becomes necessary. Nonetheless, the direct extension of IFD-BPM to the 3-D case by taking finite difference in both transverse directions squares the total number of calculations [32], leading to intolerably large computation time. Therefore, this paper proposes a scheme to handle the 3-D case without increasing calculation number. For propagation beams with circular or elliptical profiles, the 3-D case can be transferred to 2-D case as long as the singularity at the origin is properly manipulated. Because most pump beams used for parametric interaction have a circular profile [33], this scheme is practi-cally useful.

The simulation is performed in the case of the DFG in peri-odic poled lithium niobate (PPLN) for both the bulk-type and the waveguide-type materials. The extension of this scheme to other second-order nonlinear interaction is straightforward and so will not be elaborated. The comparison between the 2-D case and the 3-D case of circular beam profiles, named quasi-3-D, is shown. This paper is organized as follows. After the intro-duction, the quasi-3-D IFD-BPM is formulated and singularity at origin is solved in Section II. Then comparison between the 2-D and quasi-3-D schemes is made in Section III. In Section IV, 0733–8724/01$10.00 © 2001 IEEE

some phenomena in DFG discovered by the quasi-3-D scheme are presented. These phenomena are different from predictions by previous methods. Section V gives the conclusion.

II. FORMULATION

Propagation of light in the presence of the second-order non-linear polarization can be described by

(1a) (1b) where

and electric field and the polarization;

and linear and the second-order susceptibilities; and permittivity and permeability in vacuum. Second-order nonlinear effects are three-photon processes, so there are in general three different frequencies involved. On the other hand, the second-order nonlinear interaction is also polar-ization dependent, so the susceptibilities have in general a tensor character. For the simplicity of mathematical formulas, proper tensor elements of the linear and the second-order susceptibili-ties and the polarization of electric fields for three different fre-quencies have already been assumed.

The total electric field composed of three monochromatic waves can be expressed in the phasor notation as

c.c. (2)

where

, , angular frequencies of the three waves and for energy conservation; , , wavevectors in vacuum;

, , reference indexes.

Applying paraxial approximation, substituting (2) into (1) and equating terms with the same frequency yields

(3a)

(3b)

(3c) For 2-D scheme, the transverse Laplacian is simply described as . For the 3-D case, it becomes . If the beam is cylindrically symmetric, is -independent; then

(4)

Therefore, the two-dimensional transverse Laplacian is reduced to a one-dimensional operator. However, the first term on the right-hand side of (4) has singularity at the origin, which re-quires special treatments for the finite-difference method and will be handled later in this section.

In the finite-difference scheme, the spatial domain is divided into small regions by placing a grid over the domain. A uni-form grid is used in this work. The step sizes along the (or ) and the directions are denoted by (or ) and with and representing the indexes along these two directions, re-spectively. For example, represents the electric field at the

point (or ).

For concise expressions, the following finite difference opera-tors are defined:

(5a) (5b) (5c) (5d) (5e) where for ,

. With these definitions, the 2-D and quasi-3-D schemes could be better described.

In the IFD-BPM [29], the resulting difference equations in-volve undetermined nonlinear source terms in the next step, so iterative schemes are required to solve this problem [34]–[36]. In the beginning, one set of solutions is obtained by the RA scheme or other methods and denoted by , which is the initial guess of the electric fields for the next steps. The way of the initial guess influences the rate of convergence, but the difference is minor. An iterative algorithm described in the fol-lowing is used for simulation:

(6a)

(6c)

(or) : Quasi-3-D

where is the iteration count and is the th iteration field. In the 2-D scheme, we use expression; and in the quasi-3-D scheme, we use expression.

In the quasi-3-D IFD-BPM, there is singularity at the origin, which is not encountered in the 2-D IFD-BPM. To avoid the singularity problem for numerical calculation, Gauss’ theorem [37] is applied. From Gauss’ theorem

with being a closed loop and being the area inside. There, placement of by changes the previous formula to

Then integrating equation (3) makes the Laplacian replaced by

. With the -independence, , so the

singularity at origin is removed. Therefore, the (7) is used for the origin point

(7) where is a small circular area centered at origin and is the closed loop of .

The corresponding numerical method for the calculation at the origin is as follows. Considering the condition around the origin , a small circular area centered at

with its radius is used. The radius of the small area could be any value, but here the step size along the direction is used as the radius. Then, an “area average value” of the electric field is substituted into the area integration of (7). The “area average value” is defined in the following. Assume that

and

(8) Do the integration in (7) with a small area . Then

(9)

Similar equations can also be derived for and . These equations are used to replace (3a)–(3c). Because the di-vision by zero is now removed, the corresponding difference equations similar to (6a)–(6c) can be obtained for numerical cal-culation at the origin.

For concise expressions, a new finite-difference operator is defined

(10)

With other operators defined in (5c)–(5e), the quasi-3-D IFD around the origin can be written as

(11a)

(11b)

(11c) where is the iteration count and is the th iteration field. can be derived by various means. Equations (11a)–(11c) are the finite-difference equations used for the origin.

III. COMPARISON OFQUASI-3-DAND2-D SCHEMES In the simulation, the wavelengths and refractive indexes of the pump, the signal, and the idler waves are 0.8594 m, 1.064 m, 4.47 m, 2.171 57, 2.1575, and 2.038 14, respectively. The period of the QPM grating is 23.17 m, and the nonlinear coef-ficient is 22 pm/V. The initial power levels of the pump, the signal, and the idler waves are 0.48, 1, and 0 W, respectively. The crystal length is 2 cm with m.

(a)

(b)

(c)

Fig. 1. The diagrams of the cross-section of two waveguides. (a) The refractive-index distribution for 2-D IFD. (b) The cross section of the waveguide for 2-D IFD (the refractive index and electric field are assumed uniform along y direction). (c) The cross section of the waveguide for quasi-3-D IFD.

First, the QPM DFG in a waveguide-type PPLN is

simu-lated. The index difference is a

typ-ical value for a Ti-diffused waveguide on LiNbO . In the 2-D scheme, the width of the waveguide is assumed to be 6 m. The guided modes of the pump and the signal waves then have al-most the same full-width half-maximum (FWHM) width of 5.94 m. In the second transverse direction, the field is assumed to uniformly distribute over a range of this FWHM width. In the quasi-3-D scheme, a cylindrical waveguide with 6 m diam-eter is considered. The diagrams of the cross sections of two waveguides are shown in Fig. 1. The computation window is 100 m with the transverse grid size m. Fig. 2 shows the idler power calculated by quasi-3-D and 2-D, respec-tively, versus the propagation distance. The figure shows that the 2-D scheme predicts the growing rate of idler power only slightly larger than the 3-D scheme. At the exit facet of the 2-cm PPLN, the conversion efficiency is 22.5% and 18.2% for the 2-D scheme and the quasi-3-D scheme, respectively. Therefore, for

Fig. 2. Idler wave power versus propagation distance in the waveguide-type PPLN, calculated by 2-D and quasi-3-D schemes, respectively.

Fig. 3. Idler wave power versus propagation distance in the bulk-type PPLN, calculated by 2-D and quasi-3-D schemes, respectively.

Fig. 4. The pump wave electric field profile versus propagation distance, calculated by 2-D and quasi-3-D schemes, respectively.

the waves well confined in the waveguide, both schemes predict very similar results.

Second, the QPM-DFG in a bulk-type PPLN is simulated. The beam shapes of the pump wave and the signal wave are as-sumed to be Gaussian at the entry facet. The 25- m beam waist is located at the middle of the crystal. The computing window is 400 m. Fig. 3 shows the idler powers calculated by the 2-D scheme and the quasi-3-D scheme, respectively, versus the prop-agation distance. It can be seen that the 2-D scheme predicts a much larger growing rate of the idler power. At the exit facet, the idler power predicted by the 2-D scheme is 3.8 times the power

Fig. 5. The experimental arrangement used in [10] for comparison with our simulation.

predicted by the quasi-3-D scheme. Because the 3-D scheme takes into account the beam divergence in both transverse direc-tions, the intensity decreases rapidly, as shown in Fig. 4. This leads to reduced conversion efficiency. The 2-D scheme may also simulate the divergence in both transverse directions if the field distribution along the second transverse direction is varied according to the beam expansion. However, because the actual beam size at FWHM is still not known before it is calculated, such estimation is difficult. In comparison, the consideration of beam divergence for both transverse directions in quasi-3-D scheme is very straightforward.

The two schemes are also compared to an experiment [10], which has the experimental arrangement shown in Fig. 5. The Nd : YAG laser beam nm) is combined with the beam from a tunable external-cavity semiconductor laser using a dichroic beamsplitter. The compound external cavity of the semiconductor laser consists of a GaAlAs tapered stripe am-plifier with a 130- m output aperture and a peak gain near 855 nm, a diffraction grating for tuning, and a single stripe semicon-ductor amplifier. The laser threshold occurs at a tapered ampli-fier of A, and the output power is 820 mW at A with 0.5 W transmitted to the QPM crystal. The pump beams are focused by an cm lens, producing a 29- m FWHM beam waist ( m) at the center of a 245- m-thick 6-mm-long bulk field poled, QPM LiNbO crystal. The z-cut crystals, with a patterned electrode on the c side, were field poled with the

QPM domain period of m.

The comparison of the experiment and the simulation is shown in Fig. 6. The DFG power at m, generated by the tunable semiconductor laser at nm, is shown as a function of the pump power product . In this experiment, the crystal length is 6 mm. The simulation is thus changed for this length. Table I shows the comparison between the parameters used in the simulation and the exper-iment. The normalized nonlinear conversion efficiencies of the 2-D scheme, quasi-3-D scheme, and the experiment are 0.059, 0.0292, and 0.015%/Wcm. The normalized nonlinear conversion efficiency of the quasi-3-D scheme is closer to the experimental value than the 2-D scheme. This is reasonable because the quasi-3-D scheme is closer to the physical reality.

Fig. 6. The 4.47-m DFG idler output power versus the product of the signal power and the pump power (the pump power is kept at 0.48 W).The calculations using 2-D and quasi-3-D schemes and the experimental result are compared.

The quasi-3-D still estimates the efficiency at almost twice of the experimental value. Considering the nonideal situations in the experiment, for example, irregular variation of the periodic structure, deviation of the pump beam from the fundamental Gaussian mode, reflection at the crystal facet, and spread of the pump wavelength, the experimental conversion is inevitably less than the ideal 3-D simulation.

IV. PHENOMENADISCOVERED BYQUASI-3-D SCHEME This section shows some phenomena discovered from the quasi-3-D IFD-BPM. In the simulation, the initial beam profiles of the pump and the signal waves are Gaussian and their waist positions are assumed to be the same.

A. The Beam-Size Variation of the Idler Wave

The beam size is defined as the diameter of the cylindrical beam at which the field amplitude decays to 1 . Fig. 7(a) and (b) shows that the beam-size variation of the idler beam along the propagation distance is different from the pump and the signal beams. When the waists of the signal and the pump beams are focused at the crystal center, the beam size of the idler beam monotonically increases, growing up slowly before

TABLE I

COMPARISON OFPARAMETERSBETWEEN THESIMULATION AND THEEXPERIMENT. THEQPM PERIODSAREDIFFERENTBECAUSE OF THEDEVIATION OF

REFRACTIVEINDEX FROM THEACTUALVALUE AS ARESULT OF THEINACCURACIES IN THESELLMEIERCOEFFICIENTS ATLONGWAVELENGTH

(a)

(b)

Fig. 7. The beam sizes of the three beams for focusing position at (a) 10 000 and (b) 20 000m.

the crystal center and then quickly increasing afterwards, as shown in Fig. 7(a). The beam size of the idler wave is influenced by two factors. First, the convergence or divergence of the pump and the signal waves forces the idler wave to behave accordingly. Second, the diffraction of the idler beam

itself always results in a beam divergence. Before the crystal center, the effect of diffraction is partially canceled out by the convergence of signal and pump beams, leading to the slow increase of the beam size of the idler beam. After the crystal center, both factors cause the idler beam to diverge, so the beam size increases rapidly. If both the pump and the signal beams are focused at the end of the crystal, the convergence effect could be stronger than the diffraction effect. Then the beam size of the idler beam decreases, as shown in Fig. 7(b), but not as fast as the pump beam and the signal beam. With a proper choice of the convergent signal and pump beams, both factors could exactly cancel out one another to maintain a constant beam size of the idler beam in the crystal.

B. Beam Profiles After Propagating for a Long Distance

The beam profiles of the pump, the signal, and the idler waves at 50 000 m are shown in Fig. 8. This beam profiles are also compared to a Gaussian shape. From Fig. 8, it is clear that the beam profiles of the pump and the signal waves remain nearly Gaussian, but the generated idler wave is not Gaussian. Therefore, predictions from the plane-wave approximation or Gaussian-beam assumption [6] should fail and numerical simu-lation for beam propagation during the nonlinear conversion is necessary.

C. Influence of the Beam Size and Waist Position

In the case of DFG in bulk-type PPLN, the output idler wave power is known to be influenced by the beam size and the waist position of the pump wave and the signal waves. Fig. 9 shows the idler output power versus the propagation distance for dif-ferent beam sizes. In this calculation, the beam waists are in the middle of the 2-cm crystal. As shown in Fig. 9, the optimal beam

(a)

(b)

(c)

Fig. 8. The beam profiles at 50 000m and comparison to a Gaussian shape: (a) the pump wave, (b) the signal wave, and (c) the idler wave. [Solid line: simulated beam profile; dotted line: Gaussian shape. Both curves are almost identical for (a) and (b).]

Fig. 9. The idler output power versus the propagation distance for different waist beam sizes.

TABLE II

THEOPTIMALw , L=Z (CRYSTALLENGTH OVEROPTIMALCONFOCAL

LENGTH),ANDRANGE OF15% OFw (w THATMAKES THEOUTPUTPOWER

LESS THAN15% OF THEOPTIMALOUTPUTPOWER) VERSUSTHREE

DIFFERENTPPLN LENGTHS

Fig. 10. The normalized idler wave power versus focusing positions of the pump and the signal waves.

size at the waist is 34 m. In this situation, the crystal length is 2.18 times the confocal length ( ). For the waist beam size between 25 and 45 m, the idler output power is still more than 95% of the optimal value. This indicates that the waist beam size is not a critical issue for DFG in PPLN. The calculation had been done for several crystal lengths. The re-sults are shown in Table II, all demonstrating the insignificance of waist-size variation.

The influence of the waist position is shown in Fig. 10. It is also known that the larger the beam size is, the less significant the waist position is. However, Fig. 10 shows that the waist posi-tion is not important even for the optimal beam size. In addiposi-tion, when the beam is focused to a size as small as 20 m, the dif-ferent waist position still causes less than 15% of output power reduction, indicating that the waist position does not have sig-nificant influence as long as the beam waist is in the crystal.

V. CONCLUSION

In conclusion, a quasi-3-D IFD-BPM is developed to model second-order nonlinear interaction in both waveguide-type and bulk-type PPLNs. The quasi-3-D IFD-BPM takes into account the beam divergence in both transverse directions. In the wave-guide-type PPLN, because there is no beam divergence, the conversion efficiency calculated by the 2-D IFD-BPM and the quasi-3-D IFD-BPM is similar. For the bulk-type PPLN, the es-timated conversion efficiency by the 2-D scheme is much larger

than that predicted by the quasi-3-D scheme. The comparison of the simulation to an experiment of DFG in bulk-type PPLN clearly shows that the quasi-3-D scheme is closer to the exper-iment. The quasi-3-D IFD-BPM also reveals some novel phe-nomena in DFG, indicating that plane-wave and Gaussian-beam assumptions are not sufficient for estimating the conversion and beam propagation in second-order nonlinear devices.

REFERENCES

[1] R. L. Byer, “Quasiphase-matched nonlinear interactions and devices,”

J. Nonlinear Opt. Phys. Mater., vol. 6, pp. 549–592, 1997.

[2] L. E. Myers, R. C. Eckardt, M. M. Fejer, and R. L. Byer, “Quasiphase-matched optical parametric oscillations in bulk peri-odically poled LiNbO ,” J. Opt. Soc. Amer. B, vol. 12, pp. 2102–2116, 1995.

[3] C. Baron, H. Cheng, and M. C. Gupta, “Domain inversion in LiTaO and LiNbO by electric field application on chemically patterned crystals,”

Appl. Phys. Lett., vol. 68, pp. 481–483, 1996.

[4] M. L. Bortz and M. M. Fejer, “Annealed proton exchanged LiNbO waveguides,” Opt. Lett., vol. 6, pp. 1844–1846, 1991.

[5] D. Eger, M. Oron, M. Kartz, and A. Zussman, “Highly efficient blue light generation in KTiOPO waveguides,” Appl. Phys. Lett., vol. 64, pp. 3208–3209, 1994.

[6] M. M. Fejer, G. A. Magel, D. H. Jundt, and R. L. Byer, “Quasiphase-matched second harmonic generation: Tuning and tolerances,” IEEE J.

Quantum Electron., vol. 28, pp. 2631–2654, 1992.

[7] J. Webjorn, S. Siala, D. K. Nam, R. G. Waarts, and R. J. Lang, “Visible laser sources based on frequency doubling in nonlinear waveguides,”

IEEE J. Quantum Electron., vol. 33, pp. 1673–1685, 1997.

[8] O. Pfister, J. S. Wells, L. Hollberh, L. Zink, D. A. Van Baak, M. D. Levenson, and W. R. Bosnberg, “Continuous-wave frequency tripling and quadrupling by simultaneous tree-wave mixing in periodically poled crystal: Application to a two-step 1.19–10.71m frequency bridge,” Opt. Lett., vol. 22, pp. 1211–1213, 1997.

[9] L. Goldberg, W. K. Burns, and R. W. McElhanon, “Wide accceptance bandwidth difference frequency generation in quasiphase-matched LiNbO ,” Appl. Phys. Lett., vol. 67, no. 20, pp. 2910–2912, 1995. [10] S. Sanders, R. J. Lang, L. E. Myers, M. M. Fejer, and R. L. Byer,

“Broadly tunable mid-IR radiation source based on difference fre-quency mixing of high power wavelength-tunable laser diodes in bulk periodically poled LiNbO ,” Electron. Lett., vol. 32, pp. 218–219, 1996.

[11] C. Q. Xu, H. Okayama, and T. Kamijoh, “Broadband multichannel wavelength conversions for optical communication systems using quasiphase matched difference frequency generation,” J. Appl. Phys., vol. 34, pp. L1543–L1545, 1995.

[12] M. L. Sundheimer, Ch. Bosshard, E. W. Van Stryland, and G. I. Stegeman, “Large nonlinear phase modulation in quasiphase-matched KTP waveguides as a result of cascaded second-order process,” Opt.

Lett., vol. 18, pp. 1397–1399, 1993.

[13] C. N. Ironside, J. S. Aitchison, and J. M. Arnold, “An all-optical switch employing the cascaded second-order nonlinear effect,” IEEE J.

Quantum Electron., vol. 29, pp. 2650–2654, 1993.

[14] M. A. Arbore, O. Marco, and M. M. Fejer, “Pulse compression during second-harmonic generation in aperiodic quasiphase-matching gratings,” Opt. Lett., vol. 22, pp. 865–867, 1997.

[15] J. D. McMullen, “Optical parametric interactions in isotropic materials using a phase-corrected stack of nonlinear dielectric plates,” J. Appl.

Phys., vol. 46, pp. 3076–3081, 1975.

[16] C. Q. Xu, H. Okayama, and M. Kawahara, “Optical frequency con-versions in nonlinear medium with periodically medulated linear and nonlinear optical parameters,” IEEE J. Quantum Electron., vol. 31, pp. 981–987, 1995.

[17] V. Mahalakshmi, M. R. Shenoy, and K. Thyagarajan, “Evolution of the intensity profile of Cerenkov second-harmonic radiation with propaga-tion distance in planar waveguides,” IEEE J. Quantum Electron., vol. 32, pp. 137–144, 1996.

[18] M. Vaya, K. Thyagarajan, and A. Kumar, “Leaky-waveguide configura-tion for quasiphase-matched second-harmonic generaconfigura-tion,” J. Opt. Soc.

Amer. B, vol. 15, pp. 1322–1328, 1998.

[19] H. J. W. M. Hoekstra, “On beam propagation methods for modeling in integrated optics,” Opt. Quantum Electron., vol. 29, pp. 157–171, 1997. [20] B. Hermansson, D. Yevick, and L. Thylen, “A propagation beam method analysis of nonlinear effects in optical waveguides,” Opt. Quantum

Elec-tron., vol. 16, pp. 525–534, 1984.

[21] K. Hayata and M. Koshiba, “Numerical study of guide-wave sum-frequency generation through second-order nonlinear parametric processes,” J. Opt. Soc. Amer. B, vol. 8, pp. 449–458, 1991.

[22] , “Three-dimentional simulation of guided-wave second-harmonic generation in the form of coherent Cerenkov radiation,” Opt. Lett., vol. 16, pp. 1563–1565, 1991.

[23] F. A. Katsriku, B. B. A. Rahman, and K. T. V. Grattan, “Numerical modeling of second harmonic generation in optical waveguides using the finite element method,” IEEE J. Quantum Electron., vol. 33, pp. 1727–1733, 1997.

[24] P. S. Weitzman and U. Osterberg, “A modified beam propagation method to model second harmonic generation in optical fibers,” IEEE

J. Quantum Electron., vol. 29, pp. 1437–1443, 1993.

[25] H. M. Masoudi and J. M. Arnold, “Modeling second-order nonlinear effects in optical waveguides using a parallel-processing beam propa-gation method,” IEEE J. Quantum Electron., vol. 31, pp. 2107–2113, 1995.

[26] , “Parallel beam propagation methed for the analysis of second har-monic generation,” IEEE Photon. Technol. Lett., vol. 7, pp. 400–402, 1995.

[27] G. J. M. Krijnen, W. Torruellas, G. I. Stegeman, H. J. W. M. Heoskstra, and P. V. Lambeck, “Optimization of second harmonic generation and nonlinear phase-shifts in the Cerenkov regime,” IEEE. J. Quantum

Elec-tron., vol. 32, pp. 729–738, 1996.

[28] H. J. W. M. Heoskstra, O. Noordman, G. J. M. Krijnen, R. K. Varshney, and E. Henselmans, “Beam-propagation method for second-harmonic generation in waveguides with birefringent materials,” J. Opt. Soc Amer.

B, vol. 14, pp. 1823–1830, 1997.

[29] H.-F. Chou, C.-F. Lin, and G.-C. Wang, “An iterative finite difference beam propagation method for modeling second-order nonlinear effects in optical waveguides,” J. Lightwave Technol., vol. 16, pp. 1686–1693, 1998.

[30] H.-F. Chou, C.-F. Lin, and S. Mou, “Comparisons of finite difference beam propagarion methods for modeling second-order nonlinear effects,” J. Lightwave Technol., vol. 17, pp. 1481–1486, 1998. [31] Y. Chung and N. Dagli, “An assessment of finite difference beam

prop-agation method,” IEEE J. Quantum Electron., vol. 26, pp. 1335–1339, 1990.

[32] H.-F. Chou, “A study of quasiphase-matched second-order nonlinear ef-fects,” Masters thesis, National Taiwan University, 1998.

[33] W. Koechner, Solid-State Laser Engineering. New York: Springer-Verlag, 1976.

[34] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,

Nu-merical Recipes in C: The Art of Scientific Computing, U.K.: Cambridge

Univ. Press, 1992.

[35] J. W. Thomas, Numerical Partial Differential Equations: Finete

Differ-ence Methods. New York: Springer-Verlag, 1995.

[36] S. Nakamura, Applied Numerical Methods with Software. Englewood Cliffs, NJ: Prentice-Hall, 1991.

[37] D. G. Zill and M. R. Cullen, Advanced Engineering Mathematics: PWS, 1992.

Shing Mou, photograph and biography not available at the time of publication.

Ching-Fuh Lin (S’89–M’92–SM’00), photograph and biography not available at the time of publication.

Hsu-Feng Chou, photograph and biography not available at the time of publi-cation.