## Abstract

The second-order nonlinear frequency conversion in three-dimensional nonlinear photonic crystals is theoretically studied using coupled wave equations. A universal theoretical model is obtained, with a unified expression combining birefringence phase match, quasi-phase match, nonlinear Raman-Nath diffraction, nonlinear Čerenkov radiation and nonlinear Bragg diffraction. They are demonstrated in the numerical simulation. With the phase-matching conditions in lower dimensions extended to three dimensions, more various phenomena can be seen and corresponding mechanisms can be explained. This research enables the control of second-harmonic generation more efficiently and has potential applications in more complicated nonlinear photonic crystals.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Second-order nonlinear frequency conversion [1,2] in nonlinear photonic crystals (NPCs) [3,4] is an important research branch of nonlinear optics, which is of great practical importance to create coherent light sources at new frequency bands. It is governed by the phase-matching requirements due to the reciprocal lattice of these crystals. Therefore, the modulation of the nonlinear susceptibility enables to engineer the spatial and spectral response in these processes. As for *χ*^{(2)} processes, there are various phenomena of second harmonic wave (SH) generated in one-dimensional (1D) and two-dimensional (2D) NPCs by different types of phase match, such as birefringence phase match (BPM) [5], quasi-phase match (QPM) [6, 7], non-colinear phase match for nonlinear Raman-Nath diffraction (NRND) [8–10], nonlinear Čerenkov radiation (NCR) [11–14], as well as nonlinear Bragg diffraction (NBD) which is a degeneration of NRND and NCR [15,16]. Besides the forward emitting phenomena, there are also backward emitting cases with low conversion efficiency like backward QPM [17–19], which might have been ignored in common cases compared with forward ones. They are simple in structure and have been well researched. If the structures become more complicated with 1D and 2D NPCs extended to three-dimensional (3D) NPCs, SH that satisfying different phase-matching conditions shall emit simultaneously and mix together.

However, the research of *χ*^{(2)} processes in 3D NPCs is merely reported both in theory and experiment due to the limitation of existing manufacture techniques [20]. Thoroughly theoretical analysis helps to discriminate different types of phase match. Moreover, the interaction of different types of phase match also needs to be further studied.

In this paper, we theoretically analysed *χ*^{(2)} processes in 3D NPCs. Starting from the coupled wave equations, a universal theoretical model with a unified expression is derived. The 1D and 2D cases are also included in the universal model, while the phase-matching conditions of NRND, NCR and NBD are extended to 3D. Hence, different types of phase match in *χ*^{(2)} processes are unified into the universal model.

## 2. Theoretical model and analysis

In 1D and 2D NPCs, *χ*^{(2)} processes in different types of phase match can be classified into five categories (Fig. 1) [5–16]. Actually, these processes can also take place in 3D NPCs, furthermore *χ*^{(2)} processes shall be more various due to an extra lattice dimensionality. For example, the NCR and NRND turn into multiple-order forms compared with 1D and 2D cases.

Starting with the derivation of theoretical model to analyse the specific situation of *χ*^{(2)} processes in 3D NPCs, we only consider the interaction between fundamental wave (FW) and SH. We assume that FW incidents into the NPC along *z*-axis and generates the SH, which can be written as *E⃗*_{1,2}(*r⃗*, *t*) = *A*_{1,2}(*r⃗*)*e*^{ik⃗1,2r⃗}*e*^{−iω1,2t}, where *A*(*r⃗*) is the complex amplitude, *k⃗* is the wave vector, *ω* is the angular frequency satisfying *ω*_{2} = 2*ω*_{1}, subscripts 1 and 2 denote FW and SH, respectively. *A*_{1}(*r⃗*) is regarded as a constant under small signal approximation condition [21]. Then the paraxial coupled wave equations can be expressed as

*g*(

*x*,

*y*,

*z*) is the modulation function describing the variation of effective second-order nonlinear coefficient

*d*

_{eff}, which is periodically changed in space.

*e*

^{−(x2+y2)/w2}is the Gaussian distribution function of FW and

*w*is the beam waist radius [22]. Taking slowly varying amplitude approximation into consideration [21] and applying Fourier transformation to Eq. (1), the generated SH can be deduced as

*I*=

*A*

^{*}

*A*denotes the intensity,

*k*and

_{x}*k*are the

_{y}*x*-component and

*y*-component of

*k*

_{2}, ${\beta}_{2}={\omega}_{2}^{2}{d}_{\text{eff}}/{k}_{2}{c}^{2}$ with the light speed

*c*in vacuum. Since the derivation is conducted under paraxial approximation, which means that ${k}_{x}^{2}+{k}_{y}^{2}\ll {k}_{2}^{2}$ with the SH exiting at a small angle (within 10 degrees) to

*z*-axis, the

*z*-component of

*k*

_{2}can expand into Taylor series as ${k}_{z}={\left[{k}_{2}^{2}-\left({k}_{x}^{2}+{k}_{y}^{2}\right)\right]}^{1/2}={k}_{2}-\left({k}_{x}^{2}+{k}_{y}^{2}\right)/2{k}_{2}-{\left({k}_{x}^{2}+{k}_{y}^{2}\right)}^{2}/8{k}_{2}^{2}-\cdots $. Apparently, ${k}_{2}-\left({k}_{x}^{2}+{k}_{y}^{2}\right)/2{k}_{2}$ is the paraxial approximation of

*k*by omitting high-order terms. To eliminate errors in large angle cases,

_{z}*k*should replace ${k}_{2}-\left({k}_{x}^{2}+{k}_{y}^{2}\right)/2{k}_{2}$ in Eq. (2).

_{z}Equation (2) describes *χ*^{(2)} processes in different kinds of NPCs. In general, the *χ*^{(2)} distribution of NPCs is periodic and independent mutually in each direction. So that *g*(*x*, *y*, *z*) can be expanded into Fourier series. Then Eq. (2) can be transformed into

*n*,

*m*,

*l*denote the orders of the Fourier series in the

*x*,

*y*,

*z*directions, respectively. Now the unified expression of

*χ*

^{(2)}processes in 3D NPCs is derived as Eq. (3), containing different types of phase match. For a certain value of

*z*in Eq. (3), it is the absolute-value term that determines the intensity of

*I*

_{2}. Factors in this term are in the form of sinc or exponential function.

*I*

_{2}reaches the maximum when these factors’ values approach 1. Considering

*G*

_{0i}= 2

*π*/Λ

*, (*

_{i}*i*=

*x*,

*y*,

*z*), is the unit reciprocal vector in

*i*direction,

*I*

_{2}shall get the maximum under the transverse conditions

*k*=

_{x}*nG*

_{0x},

*k*=

_{y}*mG*

_{0y}and the longitudinal conditions

*k*= 2

_{z}*k*

_{1}−

*lG*

_{0z}, which are similar with the phase-matching conditions in 1D and 2D cases.

In reciprocal space as shown in Fig. 2(a), the phase-matching conditions in 3D NPCs can be written as ${\overrightarrow{k}}_{2}=\left(2{\overrightarrow{k}}_{1}+{\overrightarrow{G}}_{z}^{(l)}\right)+{\overrightarrow{G}}_{r}^{(n,m)}$, where ${\overrightarrow{G}}_{r}^{(n,m)}={\overrightarrow{G}}_{x}^{(n)}+{\overrightarrow{G}}_{y}^{(m)}$ [Fig. 2(b)] and ${\overrightarrow{G}}_{z}^{(l)}$ are the reciprocal vectors in transverse and longitudinal directions, respectively. Considering in each orthogonal direction, actually the transverse condition *k*_{2} sin *θ* = *G _{r}* corresponds to NRND as shown in Fig 2(c). While the longitudinal condition

*k*

_{2}cos

*θ*= 2

*k*

_{1}−

*G*corresponds to NCR as shown in Fig. 2(d). When they are both satisfied, NBD is generated. Moreover, different

_{z}*n*’s,

*m*’s,

*l*’s correspond to different orders, which means that the generated SH includes multiple orders in 3D NPCs. Especially for multiple-order NCRs, there is an extra compensation of

*G*compared with the conventional NCR having no compensation of reciprocal vectors. In addition, even though

_{z}*n*’s,

*m*’s,

*l*’s range from (−∞, ∞), the emitted SH is still of limited orders. Because there is a coefficient

*C*in Eq. (3) that determined by the specific NPC’s

_{nml}*χ*

^{(2)}structure. There are generally only several lower orders that can emit.

In NCR processes, the angle *θ* between *k⃗*_{2} and *z*-axis is the nonlinear Čerenkov angle satisfying cos*θ* = *k _{z}*/

*k*

_{2}= (2

*k*

_{1}−

*G*)/

_{z}*k*

_{2}. When Čerenkov angle goes to zero in special cases, or NBD only occurs in 1D case, it just reduces to 1D QPM form. Similarly, BPM processes could also be considered as a special case of QPM or NBD. With the addition of multiple-order SH in 3D case, different types of phase match in

*χ*

^{(2)}processes have united into a unified phase-matching condition included in Eq. (3).

## 3. Numerical simulation and discussion

Lithium niobate crystals are taken as sample to simulate the *χ*^{(2)} processes in typical 3D NPCs. The wavelengths of FW and SH are *λ*_{1} = 0.8 *μ*m and *λ*_{2} = 0.4 *μ*m, respectively. Beam waist radius of the input Gaussian beam is *w* = 50 *μ*m. And the interaction distance is *z* = 1000 *μ*m. The corresponding refractive indexes can be calculated using Sellmeier equation [23]. These parameters stay the same in all the following simulations unless otherwise stated. Two typical 3D NPCs are shown in Figs. 3(a) and 3(c). The dark areas represent positively poled domains where the modulation function *g*(*x*, *y*, *z*) = 1, while the light areas represent negatively poled domains where *g*(*x*, *y*, *z*) = −1. The duty cycle is 1:1. Hence, *C _{nml}* =

*C*with

_{n}C_{m}C_{l}*C*= 2/

_{i}*iπ*· sin(

*iπ*/2), (

*i*=

*n*,

*m*,

*l*). Between two adjacent odd orders, coefficient of the higher order is much smaller than that of the lower one. It approaches 0 with the order increasing. Besides,

*C*

_{0}= 1 and coefficients of non-zero even orders always equal to 0. As a result, 0-order is always the maximum one over other orders and only several lower positive and negative odd orders remain.

In Fig. 3(a), a polar coordinate is used for 3D cylindrical structure [24]. Specially, Λ* _{x}* and Λ

*degenerate to Λ*

_{y}*and ${\overrightarrow{G}}_{r}^{(n,m)}$degenerates to ${\overrightarrow{G}}_{r}^{({n}^{\prime})}$. The wave vector along*

_{r}*r*direction satisfies

*k*

_{r}^{2}=

*k*

_{x}^{2}+

*k*

_{y}^{2}. Λ

*is calculated as Λ*

_{z}*= 2*

_{z}*π*/(

*k*

_{2}− 2

*k*

_{1}) = 5.28

*μ*m to obtain

*I*

_{2}continuous increasement along

*z*-axis. As the

*I*

_{2}simulation shown in Fig. 3(b), the multiple-order SH displays as rings. For NRNDs, the 0-order is in the center with the positive and negative higher orders on the outer side. While for NCRs, the −1st-order, 0-order, 1st-order, 3rd-order, and higher positive orders are arranged in sequence from inside to outside. Non-zero even orders and the negative orders higher than −1 do not exist. Apparently the intensities of lower orders are stronger than those of higher orders, and the intensities of NCR rings are stronger than those of NRND nearby. Besides, the distances between adjacent orders become shorter with Λ increasing. And the position of the brightest 0-order NCR ring, which is actually the conventional NCR, is not affected by Λ

*.*

_{z}In Fig. 3(c), a Cartesian coordinate is used for cubical 3D NPC, and Fig. 3(d) is the corresponding *I*_{2} simulation result. The multiple-order NRNDs display as square spot arrays, while the multiple-order NCRs still display as rings but consist of discrete spots. Because the transverse conditions and longitudinal conditions dominate *I*_{2} simultaneously in Eq. (3), the intensities of NRNDs and NCRs are mutually interactive and restrictive. When the transverse and longitudinal conditions are both satisfied, NBDs are generated as brighter spots in NCR rings.

Since NBD is a degenerated phenomenon of NRND and NCR, it is able to get special order NCR’s enhancement by modifying the lattice period. Taking cylindrical structure for example, we keep Λ* _{z}* = 5.28

*μ*m unchanged and adjust the value of Λ

*based on the unified phase-matching condition. As shown in Fig. 4, the 0-order, 1st-order, 3rd-order NCR rings get enhanced successively with Λ*

_{r}*decreasing. In Fig. 4(a), NBD is generated at the position of 0-order NCR ring, which becomes much brighter as the strongest one over other orders. Keeping Λ*

_{r}*decreasing, NBD is generated at the position of 1st-order NCR, and it is almost enhanced to the level of the adjacent lower order as shown in Fig. 4(b). So as the 3rd-order NCR enhancement shown in Fig. 4(c).*

_{r}## 4. Dimension reduction analysis in the universal model

Different (*n*, *m*, *l*)’s correspond to different dimensional NPCs. They are all contained in Eq. (3). Hence, 1D and 2D cases can also be analysed.

When (*n*, *m*, *l*) = (0, 0, 0) describing *χ*^{(2)} processes in bulk crystals with no periodic structure, Eq. (3) reduces into

When (*n*, *m*, *l*) = (0, 0, *l*), it means there is no periodic structure in *x* or *y* direction. The standard 3D NPC [Fig. 5(a)] reduces to 1D case [Fig. 5(b)] with Eq. (3) reducing into

*χ*

^{(2)}processes in 1D NPCs. The simulation results in Figs. 5(c) and 5(d) exactly show the conventional QPM processes [7]. Besides, Eq. (5) also gives the distribution of generated SH in transverse directions. If the incident direction changes as shown in Fig. 5(e), we get (

*n*,

*m*,

*l*) = (

*n*, 0, 0). Then Eq. (3) reduces into

*x*-direction as shown in Fig. 5(f).

When (*n*, *m*, *l*) = (*n*, *m*, 0), it means there is no periodic structure in *z* direction. The 3D NPC reduces to 2D case [Fig. 6(a)] with Eq. (3) reducing into

*g*(

*x*,

*y*,

*z*).

Similarly, we can obtain the SH intensity distribution in more *χ*^{(2)} structures, and patterns can be designed to realize desired SH distribution. Furthermore, making appropriate artificial *χ*^{(2)} structures in NPCs would allow us to control the behavior of harmonic generation more efficiently, which strongly extends potential applications in optical control, harmonic generation, as well as in photonics studies.

## 5. Conclusion

In this paper, a theoretical study of *χ*^{(2)} processes in 3D NPCs is taken, and a universal model with a unified expression that combines birefringence phase match, quasi-phase match, nonlinear Raman-Nath diffraction, nonlinear Čerenkov radiation and nonlinear Bragg diffraction in 1D, 2D and 3D NPCs is obtained. The phase-matching conditions in lower dimensions are extended to three dimensions, and multiple-order nonlinear Raman-Nath diffractions and nonlinear Čerenkov radiations are generated. The numerical simulation of 3D cases is demonstrated, as well as 1D and 2D cases are re-deduced with the universal model. They agree well with previous researches. Furthermore, other structures of NPCs can be designed to generate desired SH, which enables the control of harmonic generation more efficiently.

## Funding

National Key R&D Program of China (2017YFA0303700); National Natural Science Foundation of China (NSFC) (11734011, 11604318); The Foundation for Development of Science and Technology of Shanghai (17JC1400400)

## References and links

**1. **M. M. Fejer, “Nonlinear optical frequency conversion,” Phys. Today **47**, 25–33 (1994). [CrossRef]

**2. **R. Lifshitz, A. Arie, and A. Bahabad, “Photonic quasicrystals for nonlinear optical frequency conversion,” Phys. Rev. Lett. **95**, 133901 (2005). [CrossRef] [PubMed]

**3. **C. M. Bowden and A. M. Zheltikov, “Nonlinear optics of photonic crystals introduction,” J. Opt. Soc. Am. B **19**, 2046–2048 (2002). [CrossRef]

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

**5. **W. Gandrud, “On the possibility of phase-matched infrared mixing using induced birefringence,” IEEE J. Quantum Electron. **7**, 132–133 (1971). [CrossRef]

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

**7. **L. E. Myers, R. Eckardt, M. Fejer, R. Byer, W. Bosenberg, and J. Pierce, “Quasi-phase-matched optical parametric oscillators in bulk periodically poled LiNbO_{3},” J. Opt. Soc. Am. B **12**, 2102–2116 (1995). [CrossRef]

**8. **S. M. Saltiel, D. N. Neshev, W. Krolikowski, A. Arie, O. Bang, and Y. S. Kivshar, “Multiorder nonlinear diffraction in frequency doubling processes,” Opt. Lett. **34**, 848–850 (2009). [CrossRef] [PubMed]

**9. **Y. Sheng, Q. Kong, W. Wang, K. Kalinowski, and W. Krolikowski, “Theoretical investigations of nonlinear Raman–Nath diffraction in the frequency doubling process,” J. Phys. B: At. Mol. Opt. Phys. **45**, 055401 (2012). [CrossRef]

**10. **A. Vyunishev, V. Slabko, I. Baturin, A. Akhmatkhanov, and V. Y. Shur, “Nonlinear Raman–Nath diffraction of femtosecond laser pulses,” Opt. Lett. **39**, 4231–4234 (2014). [CrossRef] [PubMed]

**11. **P. Tien, R. Ulrich, and R. Martin, “Optical second harmonic generation in form of coherent Cerenkov radiation from a thin-film waveguide,” Appl. Phys. Lett. **17**, 447–450 (1970). [CrossRef]

**12. **Y. Sheng, Q. Kong, V. Roppo, K. Kalinowski, Q. Wang, C. Cojocaru, and W. Krolikowski, “Theoretical study of Čerenkov-type second-harmonic generation in periodically poled ferroelectric crystals,” J. Opt. Soc. Am. B **29**, 312–318 (2012). [CrossRef]

**13. **H. Ren, X. Deng, Y. Zheng, N. An, and X. Chen, “Nonlinear Cherenkov radiation in an anomalous dispersive medium,” Phys. Rev. Lett. **108**, 223901 (2012). [CrossRef] [PubMed]

**14. **X. Zhao, Y. Zheng, H. Ren, N. An, X. Deng, and X. Chen, “Nonlinear Cherenkov radiation at the interface of two different nonlinear media,” Opt. Express **24**, 12825–12830 (2016). [CrossRef] [PubMed]

**15. **G. Pan, R. Kesavamoorthy, and S. A. Asher, “Optically nonlinear Bragg diffracting nanosecond optical switches,” Phys. Rev. Lett. **78**, 3860 (1997). [CrossRef]

**16. **S. M. Saltiel, D. N. Neshev, R. Fischer, W. Krolikowski, A. Arie, and Y. S. Kivshar, “Generation of second-harmonic conical waves via nonlinear Bragg diffraction,” Phys. Rev. Lett. **100**, 103902 (2008). [CrossRef] [PubMed]

**17. **A. C. Busacca, S. Stivala, L. Curcio, A. Tomasino, and G. Assanto, “Backward frequency doubling of near infrared picosecond pulses,” Opt. Express **22**, 7544–7549 (2014). [CrossRef] [PubMed]

**18. **Y. J. Ding, J. U. Kang, and J. B. Khurgin, “Theory of backward second-harmonic and third-harmonic generation using laser pulses in quasi-phase-matched second-order nonlinear medium,” IEEE J. Quant. Electron. **34**, 966–974 (1998). [CrossRef]

**19. **M. Lauritano, A. Parini, G. Bellanca, S. Trillo, M. Conforti, A. Locatelli, and C. D. Angelis, “Bistability, limiting, and self-pulsing in backward second-harmonic generation: a time-domain approach,” J. Opt. A: Pure Appl. Opt. **8**, S494–S501 (2006). [CrossRef]

**20. **P. Ferraro, S. Grilli, and P. D. Natale, *Ferroelectric Crystals for Photonic Applications: Including Nanoscale Fabrication and Characterization Techniques*, vol. 91 (Springer Science & Business Media, 2013).

**21. **R. W. Boyd, *Nonlinear Optics* (Academic, 2003).

**22. **J. Wen and M. Breazeale, “A diffraction beam field expressed as the superposition of Gaussian beams,” J. Acoust. Soc. Am. **83**, 1752–1756 (1988). [CrossRef]

**23. **O. Gayer, Z. Sacks, E. Galun, and A. Arie, “Temperature and wavelength dependent refractive index equations for MgO-doped congruent and stoichiometric LiNbO_{3},” Appl. Phys. B **91**, 343–348 (2008). [CrossRef]

**24. **J. Chen and X. Chen, “Generation of conical and spherical second harmonics in three-dimensional nonlinear photonic crystals with radial symmetry,” J. Opt. Soc. Am. B **28**, 241–246 (2011). [CrossRef]

**25. **A. Arie, N. Habshoosh, and A. Bahabad, “Quasi phase matching in two-dimensional nonlinear photonic crystals,” Opt. Quant. Electron. **39**, 361–375 (2007). [CrossRef]

**26. **X. Wang, X. Zhao, Y. Zheng, and X. Chen, “Theoretical study on second-harmonic generation in two-dimensional nonlinear photonic crystals,” Appl. Opt. **56**, 750–754 (2017). [CrossRef] [PubMed]

**27. **S. Inoue and Y. Aoyagi, “Design and fabrication of two-dimensional photonic crystals with predetermined nonlinear optical properties,” Phys. Rev. Lett. **94**, 103904 (2005). [CrossRef] [PubMed]

**28. **A. Arie and N. Voloch, “Periodic, quasi-periodic, and random quadratic nonlinear photonic crystals,” Laser Photonic Rev. **4**, 355–373 (2010). [CrossRef]

**29. **N. Broderick, G. Ross, H. Offerhaus, D. Richardson, and D. Hanna, “Hexagonally poled lithium niobate: a two-dimensional nonlinear photonic crystal,” Phys. Rev. Lett. **84**, 4345 (2000). [CrossRef] [PubMed]