Numerical investigation of effective nonlinear coefficient model for coupled third harmonic generation

In this paper, the optimal solution of effective nonlinear coefficient of quasi-phase-matching (QPM) crystals for coupled third harmonic generation (CTHG) was numerically investigated. The effective nonlinear coefficient of CTHG was converted to an Ising model for optimizing domain length distributions of aperiodically poled lithium niobate (APPLN) crystals with lengths as 0.5 mm and 1 mm, and fundamental wavelengths ranging from 1000 nm to 6000 nm. A method for reconstructing crystal domain poling weight curve of coupled nonlinear processes was also proposed, which demonstrated the optimal conversion ratio between two coupled nonlinear processes at each place along the crystal. In addition, by applying the semidefinite programming, the upper bound on the effective nonlinear coefficients deff for different fundamental wavelengths were calculated. The research can be extended to any coupled dual \c{hi}(2) process and will help us to understand better the dynamics of coupled nonlinear interactions based on QPM crystals.


Introduction
As one of the most important inventions of our time, laser technology has very important applications in numerous fields, such as precision measurement [1,2], laser processing [3], fundamental physics [4,5], and quantum information [6,7].However, the wavelengths of lasers depend on the energy level structure of gain mediums, which occupy limited ranges with respect to the whole optical spectrum.Nonlinear optical frequency conversion techniques are essential to extend the applications [8].Quasi-phase matching (QPM) [9][10][11] is a widely used conversion technique, in which, the sign of the second-order susceptibility is conventionally changed periodically in a piece of nonlinear crystal (lithium niobate, for example) to obtain a periodically poled structure, offering an additional reciprocal lattice vector to satisfy the momentum conservation condition.Generally, a periodically poled lithium niobate (PPLN) crystal can only satisfy a single χ (2) process in a narrow bandwidth.However, in practical applications, multi-process frequency conversion in a single QPM crystal is expected for rich spectrum availability with high efficiency, low cost and decent stability, which cannot be realized by conventional PPLN crystals.
To solve this problem, successive researchers have broken through the limitation of periodically poling design and proposed quasi periodically poled lithium niobate (QPPLN) [12], chirped periodically poled lithium niobate (CPPLN) [13,14], and aperiodically poled lithium niobate (APPLN) [15,16] crystals.Through flexible modulation of the domain lengths, multiple nonlinear processes and wide spectrum bandwidth can be realized in a single QPM crystal.For multi-process frequency conversion, APPLN can achieve higher conversion efficiency than conventional cascaded PPLNs [26].Optimizing the poling direction and the length of the each domain in APPLN can be obtained by maximizing the effective nonlinear coefficient deff through methods such as simulated annealing (SA), genetic algorithm [17] and even quantum annealing [18].Taking the CTHG [9,15] that contains both second harmonic generation (SHG, ω+ω=2ω) as well as sum frequency generation (SFG,ω+2ω=3ω) as an example, for different fundamental wavelengths λ1, what are the distribution modes of the poling directions of the corresponding optimal CTHG crystal?How does the effective nonlinear coefficient vary with the fundamental wavelengths?Does deff have a theoretical upper bound?These have been less addressed in previous studies.The understanding of these questions will help us to understand the conversion mechanism of CTHG better.
In this paper, we convert the effective nonlinear coefficient model to the Ising model [19], and use SA algorithm to search the optimal domain distributions of APPLN and calculate the corresponding CTHG deff with crystal lengths as 0.5mm and 1mm , respectively, when fundamental wavelength varys from 1000nm to 6000nm.We also propose a method to reconstruct the poling weight curves of domain distributions from the obtained APPLNs.In addition, by applying the semidefinite programming [21], we calculated the upper bound of the deff for different fundamental wavelengths.

Methodology for CTHG crystal design
The QPM technique introduces reciprocal lattice vectors by changing the sign of the crystal's nonlinear coefficients within a coherence length, which compensates the phase mismatch during nonlinear frequency conversion process.As shown in Fig. 1, the APPLN crystal is divided into equal-sized domains, each can have either "up" or "down" poling direction.(It should be noted that in some articles, the term "domain" have different definition which refers to the length of continuous poling in the same direction, in that case each domain may be unequal-sized.)Optimized arrangement of the directions of the domains can achieve the most efficient nonlinear frequency conversion.The CTHG process contains both SHG and SFG nonlinear processes.Under plane-wave approximation, the coupled wave equations [22] are where i is imaginary unit, E1, E2 and E3 are complex amplitudes of fundamental, second and third harmony, respectively.Δk1=n2ω2/c-2n1ω1/c and Δk2=n3ω3/c-n1ω1/c-n2ω2/c are the phase mismatches of the SHG and SFG processes, respectively.d33 is the nonlinear coefficient of the lithium niobate crystal in the 33 direction, c is the speed of light in a vacuum, and ω1, ω2, and ω3 are the frequencies of the fundamental wave, second harmonic wave, and third harmonic wave, respectively.n1, n2, and n3 are the corresponding refractive indices, respectively.All simulations performed at 80℃ in this paper.d(z) is the poling direction of the domain at position z, where d(z)= +1 means "up" and -1 means "down".Under the small signal approximation, the CTHG conversion efficiency is given by [15]   where I1 is the fundamental wave intensity, I3 the third harmonic wave intensity, L the crystal length, λ the fundamental wavelength, the vacuum permittivity, and deff the effective nonlinear coefficient defined as [15] 2 1 2 0 0 2 ( ) ( ) For given pump and signal wavelengths, Δk1 and Δk2 can be obtained from Sellmeier's equation [23].Equation (2) shows that to maximize ηTHG, we need to find an optimal set of domain arrangements d(z) that maximizes the deff.The crystal is divided into N domain with equal length Δx = L/N.Discretize the equation (3) to obtain ( 1) Equation ( 4) can be abbreviated as Here, is the poling direction of the number m crystal domain, and J(m,n) is the coupling coefficient which is expressed as Simulations show that for different m d , complex value of   starts at zero in the complex plane, when the maximum mode lengths in different directions are approximately equal.To simplify the computation, we define the objective function f as the projection of deff on the real number axis, where the real coupling coefficient The Ising model [19,24] plays an important role in statistical physics [27], which defined as , where spins and coupling coefficient .
From equation (7) , it can be seen that the objective function f is an Ising model [19,24], whose value is determined by the poling direction of each crystal domain, and the maximum value can be obtained by using SA algorithm.
For an APPLN with N calculation crystal domains, there are 2 N different combinations, and there is no effective algorithm yet to find the global optimal solution in reasonable time using classical computers.Article [18] broke through the limitation of classical computers and proposed a method of quantum annealing algorithm to search for the global optimal solution, which is a very promising algorithm in the field of APPLN design.However, the current quantum computers have limited number of quantum bits, which hampers the accuracy of crystal domain design.For these reasons, this paper still chooses SA algorithm to search for the optimal solution.In this way, the solution might be a local optimum, but it is the most feasible and effective solution under the current computer development.Unless otherwise specified, "optimal" in this paper refers to the best result obtained by our search in a reasonable time.
For the fundamental wavelength λ1, its second harmonic wavelength and third harmonic wavelength are , respectively.The phase mismatch are for SHG and SFG, respectively, where 2 is refractive indices (i=1,2,3).

Optimal deff at different wavelengths
In our simulation, the fundamental wavelength ranges from 1000nm to 6000nm with a 10nm interval for crystal length L as 0.5mm and 1mm, respectively.L to be limited in 1mm is benefit to computing time and reasonable for the applications in femtosecond region.The SA simulated maximum value of deff is shown in Fig. 2.
Fig. 2. The values of the deff corresponding to the optimal APPLN crystal for fundamental wavelengths of 1000nm-6000 nm when the crystal length L are 0.5 mm (blue) and 1 mm (red), respectively.
In Fig. 2, the effective nonlinear coefficient deff at most of fundamental wavelengths are approximately equal to a fixed value of about 0.235, and the curves show peaks at the wavelengths near 1420nm, 2080nm, and 3550nm.The largest peak appears near λ1 = 3550nm.In order to understand the mechanism, we calculated the phase mismatches Δk1 and Δk2 for SHG and SFG at different fundamental wavelengths, and the results show that the peaks of the curves of deff in Fig. 2 are locating at the positions where  At the degenerate point the optimal APPLN is actually a PPLN crystal.We will prove this mathematically in the following.When Δk1=Δk2, the objective function ( 7) is transformed into Using the identical equation where, 10) can be written as From above equation, 1 ( , , ) Apply triangular equation cos cos 2sin sin 2 2 The equal sign of equation ( 13) holds if and only if the poling direction of the number m crystal domain is Hence the APPLN is degenerated to PPLN with a poling period of  The optimal crystal domain distributions, searched by SA algorithm, for Δk2/Δk1 = 3, 2, 1 are illustrated in Fig. 5(a), (b), (c), respectively.In Fig. 5(a), the poling length of the first part of the crystal is approximately equal to the period for SHG process, which is the case discussed in above paragraph.This is the reason why the deff corresponding to the optimal APPLN crystal near the fundamental wavelength of 1420 nm appears a peak in Fig. 2.

Optimal poling weight analysis
Different from a conventional PPLN for a single second-order nonlinear process with only one reciprocal lattice vector Δk, an APPLN for CTHG needs two reciprocal lattice vectors, Δk1 and Δk2.The efficiencies of SHG and SFG are coupled to each other, and each process happens anywhere in the crystal, so there must be an optimal balance in the superlattice design strategy for the ratio between these two processes at each place of the crystal along the propagation direction.Therefore, we define the poling function of the crystal as Thesis [20] gives formula for the optimal distribution of the weight ratio c(m) under certain approximations (see Supplementary Material).Fig. 6 shows the comparison of deff between the optimal APPLN calculated by Zhang's formula and those by the SA algorithm when the length of the crystal is L=0.5mm and 1mm.In Fig. 6, the curves by Zhang's formula are very close to our results in most positions, but have obvious deviation in the vicinities of λ1=1420nm and λ1=3550nm.In these areas the APPLN superlattice arrangements searched by SA have higher deff than those by Zhang's formula.We will discuss these through analysis and simulation in the following.The basic idea of our calculation is to extract the sum of the terms related to the crystal domain dm in the objective function (7), denoted as fm, and by calculation we show that fm can be transformed into the form of the poling function (14) .In order to keep the objective function as close to the maximum as possible, fm > 0 is needed to obtain the poling direction of dm.From the expression of dm at this point, we can derive reversely the distribution curve of the weights a(m) and b(m) corresponding to the actual optimal APPLN.
The objective function ( 7) can be rewritten as where ( ) cos( Neglecting the coefficient 2 , sum of the terms in ( 15) associated with dm is Then ( 16) can be written as which can be in the form of a cosine function If the sign of the m th domain then fm > 0, and the objective function f is closer to the maximum.Since m is an arbitrary integer between 1 ~N , equation (20) presents the solution for every domain in the crystal.
Comparison of ( 14) and (20) shows that weight factors of the poling period for SHG and SFG are and initial phases in ( 14) are In section 3, the optimal APPLN for the fundamental wavelength from 1000nm to 6000nm were calculated by SA.Applying formula (21) , we can derive reversely the distribution of weights a(m) , b(m) and c(m) from the optimal APPLN designs at λ1 = 1000nm, 3300nm, 3550nm, respectively.The results are shown in Fig. 7.The optimized domain distributions of the APPLNs calculated from Zhang's formula and SA algorithm at above three fundamental wavelengths are illustrated in Fig. 7.
As shown in Fig. 7, the weight ratio curve c(m) at λ1 = 1000nm is very smooth and close to that from Zhang's formula (a1), but they significantly deviated from each other (a2), corresponding to different domain distributions as shown in (b2), (c2), respectively.Because the domain distribution in Fig. 7(c2) was optimized by SA and the corresponding deff obtained by ( 7) is higher than that from Zhang's formula as shown in Fig. 6(a) at 3300nm.We think the APPLN designed by us is better than Zhang's formula at this fundamental wavelength.At degenerate point of 3550nm, the solid red line of weight ratio curve c(m) as shown in Fig. 7(a3) is a straight line with negative slope, which is also different from the solid blue line, but both are corresponding to almost the same domain distributions in (b3), (c3), respectively.This is because that at degenerate point, Δk1=Δk2, any weight ratio curve gives the result, a PPLN.It is easy to expect different deff curves of different nonlinear materials.We plot the simulation results in Fig. 8 for APPLN and aperiodically poled lithium tantalate (APPLT).It shows the curve of APPLT shifts to the short wavelength with degenerate point at 3370nm.
Actually, the method we proposed above is not only suitable for CTHG (SHG+SFG), but also any coupled two


processes (with corresponding Δk1 and Δk2) in a single QPM crystal, such as SHG+DFG, DFG+SFG, etc. Fig. 8. Values of the effective nonlinear coefficients corresponding to the optimal APPLN crystal (blue) and APPLT crystal (red) for fundamental wavelengths from 1000 to 6000 nm when the crystal length L is 0.5 mm.It can be seen that the effective nonlinear coefficients of the two crystal obey the same distribution law despite the different types of crystal.

Upper bond of deff
In previous sections, we used a SA algorithm to calculate the "optimal" CTHG crystal for different fundamental wavelengths, however, such a search may be limited to the local optimal solution rather than the global optimal solution.One may naturally ask that how far the local optimal solution searched by the SA method is away from the global optimal.
The direct way to solve this problem is to compute the global optimal solution of the objective function (7) .However, this is very difficult.For a crystal with N domains (calculation steps), there are 2 N different arrangements, resulting in no guaranteed global optimal searched out in a reasonable time by conventional algorithm.Therefore, we can only use an indirect method to estimate the upper bound of the optimal solution of the objective function (7).
Here, we use the semidefinite programming (SDP) [21] to compute the upper bound of the objective function f .The basic idea is to cancel the restriction of independent variable dm in equation (7) only to be +1 and -1, and "relax" dm into an N-dimensional vector that satisfies |xm|=1.By applying the semidefinite programming, the original problem can be transformed into , m . ., where X=[x1, x2, ..., xN] is a N×N semidefinite positive matrix construct by relaxation vector xm .
The problem (23) after relaxation is a polynomial time problem.Its globally optimal can be found in a reasonable time.This is the essential difference between the relaxed problem (23) and the original problem (22).Obviously, the relaxation problem (23) extends the scope of the definition domain of the original problem (22).But if the vectors in problem ( 23) are restricted to ( , 0, 0, , 0) m q    , then the relaxation problem (23) degenerates to the original problem (22).For a crystal with L=0.5 mm, N=500 when Δx=1μm.The global optimal solution of the problem (23) , which is the theoretical upper bound of the original problem (22) , can be derived (see Fig. 9).It shows that the optimal solutions derived by the SA are very close to the theoretical upper bound.It is worth noting that the red curve in Fig. 9 is different from that in Fig. 2 at short wavelengths, since here we use Δx=1μm, which is not fine enough to get high deff comparing to Δx=0.1μm used in Fig. 2, but benefit to reduce computing time.
The N-dimensional vectors corresponding to the optimal solutions derived from equation ( 23) are nearly evenly distributed on a unit circle in a two-dimensional plane when the wavelengths are far away from the degenerate point (λ1=3550), i.e. numerical result shows that each xm corresponding to the solution of problem ( 23 one of the cases, when the fundamental wavelength λ1 = 1200 nm.However, there is no such property in the vicinity of the degenerate point (Fig. 9 (c)), the reason of which is to be investigated in the future.

Conclusion
This paper focused on coupled dual second-order nonlinear process and converts the effective nonlinear coefficient into an Ising model, and used the simulated annealing method to find the domain design of the aperiodically poled QPM crystals with the largest effective nonlinear coefficient.Taking the CTHG as an example, the domain distribution of the optimal APPLN is calculated for the fundamental frequency light in the range from 1000nm to 6000nm.The simulation results show that the deff curve has peaks at Δk2/Δk1=1, 2, and 3, respectively, over the entire calculated spectral range.We demonstrate that the optimal APPLN degenerates into a PPLN at the degenerate point (Δk2=Δk1), where the deff can be as high as 0.41, which is 1.7 times higher than that at the general wavelength.In this paper, we proposed a reconstruction of the weight ratio of the dual nonlinear process at different positions within the crystal based on the optimized APPLN domain distribution.The obtained weight curves were compared with those generated from the analytical approximation formula.The numerical solutions and the approximate analytical curves are quite different on both sides of the degenerate point.We also estimated the theoretical upper bound of deff using the semidefinite programming, and the results show that the optimal solution obtained by the SA method is close to the theoretical upper bound, especially in the vicinity of the degenerate point.The methods proposed in this paper are applicable to any coupled dual second-order nonlinear process, provide guidance for efficient coupled QPM crystal design, and contributes to the understanding of the dynamics of multiprocess nonlinear frequency conversion.

S1. Histogram representation of APPLN crystal poled domains
In this paper, we used histograms to represent the domain distributions of APPLNs, which is intuitive to illustrate the change of the length of the crystal domains along the light propagation direction.
In our simulation, we used evenly spaced steps(Δx) for maximizing the effective nonlinear coefficient(deff) by determining the direction of each step.As illustrated in Fig. S1(a), a simple APPLN model contains 14 crystal segments with Δx=1μm (we used shorter steps in our real simulation).Assuming the result of optimal poling direction is as Fig. S1(a), then the corresponding domain distribution is as Fig. S1(b), and the histogram is Fig.S1(c).The red color indicates that the poled direction is upward and the blue color downward.

S2. Derivation of coupled third harmony generation small-signal solutions
The coupled third harmony generation (CTHG) process contains both SHG and SFG nonlinear processes, whose coupled wave equations are where i is imaginary unit, Δk1 and Δk2 are the phase mismatches of the SHG and SFG processes, respectively.d33 is the nonlinear coefficient of the lithium niobate crystal in the 33 direction, c is the speed of light in a vacuum, and ω1, ω2, and ω3 are the frequencies of the fundamental wave, second harmonic wave, and third harmonic wave, respectively.n1, n2, and n3 are the corresponding refractive indices, respectively.d(z) is the poling direction of the domain at position z, where d(z)= +1 means "up" and -1 means "down".In the small-signal approximation, assuming that the nonlinear conversion is not too strong and that the fundamental frequency light E1 remains constant in the crystal, |E3| is very small compared to |E2|, and the second term on the right side of (S2) can be neglected, E2 can be obtained by The effective nonlinear coefficient is expressed as

Fig. 1 .
Fig. 1.Schematic diagram of the APPLN structure.A lithium niobate crystal is divided into multiple equal-sized domains, and each domain takes either "up" or "down" poling direction.The consecutive domains with the same direction form a real poling domain for QPM crystal manufacturing process.
be understood that since Δk1=Δk2 at the degenerate point, a PPLN with identical poling periods of two processes simultaneously.For understanding the mechanism of CTHG deff enhancement at λ1 = 1420nm, where Δk2/Δk1 = 3, we visualize the local gains of λ2 and λ3 when they propagate in unpoled and poled nonlinear crystals as shown in Fig.4(a) and (b), respectively.Their accumulated gains are shown Fig.4(c) and (d), respectively.Fig.4(d)indicates that the crystal poled for SHG also present gain for SFG.In details, the accumulated gain of SHG in a single poling period 1  is continuously increasing, meanwhile that of SFG also have some increasement after one poling period 1  because of the inversion of nonlinearity after 1

Fig. 4 .
Fig.4.At the fundamental wavelength of 1420nm, Δk2/Δk1=3, when the crystal is poled with a period of 2π/Δk1, the SHG can be one-order phase matched and the SFG can be third-order phase matched.

Fig. 5 .
Fig. 5. Optimal CTHG crystal with fundamental wavelength of 1420nm(a), 2080nm(b), and 3550nm(c), with crystal length as 0.5 mm.The blue bar indicates the poling direction downward and the red bar indicates upward, as detailed in Supplementary Material.
) where m denotes the number of the crystal domain in the beam direction, a(m) and b(m) are the weight factors of the poling period for SHG and SFG, respectively.1 the proportion of SHG in the APPLN poling.

Fig. 6 .
Fig. 6.Comparison of the deff calculated from Zhang's formula(blue) and those by our SA simulation (red) with crystal length L = 0.5 mm (a) and 1 mm (b).

Fig. 9 .
Fig. 9. (a) Maximum deff obtained by the SA (red) and the theoretical upper bound derived by the SDP (blue).N vectors plotted in two-dimensional plane corresponding to the optimal solution of equation(23) at λ1 = 1200 nm (b) and λ1 = 3400 nm (c).

Funding.
National Key Research and Development Program of China (2022YFB4601103); Guangdong Basic and Applied Basic Research Foundation (2020B1515120041).Disclosures.The authors declare no conflicts of interest.

Fig. S1 .
Fig. S1.Histogram representation of the length distribution of APPLN crystal domains.(a) Optimal poled direction of each crystal domain (b) Alternative equivalent representation of the poled direction of crystal domains (c) Histogram representation of the optimal poled direction Zhang[1]  derived the poled period distribution of the CTHG optimal APPLN under certain approximations in Section 4.2 of his Ph.D. thesis (C.Zhang, "Dynamic theoretical study of the coupling phase matching process of nonlinear optical superlattices," PhD thesis of NanJing University, 2004.(written in Chinese)), and we restate his results here.Assuming that the polarization period of the crystal is varying weight function k(x) is introduced to regulate the local Fourier coefficients at x. Assume that k(x) is a slow-varying function with respect to the crystal length L. It can be shown that the local Fourier coefficients at x contribute efficiently to the overall CTHG as