1. 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–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 χ⁽²⁾ 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 d_eff 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 λ₁, 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 d_eff 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 d_eff with crystal lengths as 0.5 mm and 1 mm, respectively, when fundamental wavelength varys from 1000 nm to 6000 nm. 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 d_eff for different fundamental wavelengths.

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

 

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.

Download Full Size | PDF

The CTHG process contains both SHG and SFG nonlinear processes. Under plane-wave approximation, the coupled wave equations [22] are

For given pump and signal wavelengths, Δk₁ and Δk₂ 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 d_eff. The crystal is divided into N domain with equal length Δx = L/N. Discretize the equation (3) to obtain

Equation (4) can be abbreviated as

Here, $d_m\equiv d(m)\in {\rm \{ } + 1,-1\}$ 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 ${d_m}$, complex value of $\sum\limits_{m = 1}^N {\sum\limits_{n = 1}^{m - 1} {J(m,n){d_m}{d_n}} } $ 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 d_eff on the real number axis,

The Ising model [19,24] plays an important role in statistical physics [27], which defined as $H ={-} \sum\limits_{ij} {{J_{ij}}{\sigma _i}{\sigma _j}}$ , where spins ${\sigma _{ij}}{\sigma _i} \in \{{ + 1, - 1} \}$ and coupling coefficient ${J_{ij}} \in {{\mathbb R}}$. 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 λ₁, its second harmonic wavelength and third harmonic wavelength are ${\lambda _2} = \frac{{{\lambda _1}}}{2}$ and ${\lambda _3} = \frac{{{\lambda _1}}}{3}$, respectively. The phase mismatch $\Delta {k_1} = {k_2} - 2{k_1}$ and $\Delta {k_2} = {k_3} - {k_1} - {k_2}$ are for SHG and SFG, respectively, where ${k_i} = \frac{{2\pi }}{{{\lambda _i}}} \cdot {n_i}$, in which n_i is refractive indices (i = 1,2,3).

3. Optimal d_eff at different wavelengths

In our simulation, the fundamental wavelength ranges from 1000 nm to 6000 nm with a 10 nm interval for crystal length L as 0.5 mm and 1 mm, respectively. L to be limited in 1 mm is benefit to computing time and reasonable for the applications in femtosecond region. The SA simulated maximum value of d_eff is shown in Fig. 2.

 

Fig. 2. The values of the d_eff 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.

Download Full Size | PDF

In Fig. 2, the effective nonlinear coefficient d_eff 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 1420 nm, 2080nm, and 3550 nm. The largest peak appears near λ₁ = 3550 nm. In order to understand the mechanism, we calculated the phase mismatches Δk₁ and Δk₂ for SHG and SFG at different fundamental wavelengths, and the results show that the peaks of the curves of d_eff in Fig. 2 are locating at the positions where $\frac{{\Delta {k_2}}}{{\Delta {k_1}}}$ are integers as show in Fig. 3. Especially, at λ₁= 3550 nm where the d_eff curves have the highest peaks, Δk₁ equals to Δk₂, which referred as degenerate point. Fig. 2 and Fig. 3 show that for high efficient CTHG, it’s better to choose the fundamental wavelengths in the cases where $\frac{{\Delta {k_2}}}{{\Delta {k_1}}}$ are integers. Of course, the locations of these special wavelengths depend on different crystal materials. It can be of benefit to obtain supercontinuum generation if the fundamental wavelengths are chosen at the degenerate point of QPM crystals [25].

 

Fig. 3. The ratio of the phase mismatch (Δk₂/Δk₁) between the SFG (ω+2ω=3ω) and the SHG (ω+ω=2ω) of the PPLN crystal versus the fundamental wavelength (=2π/ω) when the temperature is 80 °C. The results show that Δk₂ is an integer multiple of Δk₁ when the fundamental wavelength λ₁= 1420 nm, 2080nm and 3550 nm, respectively.

Download Full Size | PDF

At the degenerate point the optimal APPLN is actually a PPLN crystal. We will prove this mathematically in the following. When Δk₁=Δk₂, the objective function (7) is transformed into

Using the identical equation ${\left( {\sum\limits_{i = 1}^N {{a_i}} } \right)^2} = 2\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^{i - 1} {{a_i}{a_j}} } + \sum\limits_{i = 1}^N {a_i^2} $, it is possible to transform $f({d_1},d{}_2, \cdot{\cdot} \cdot ,{d_N})$ into

From above equation, $f({d_1}, \cdot{\cdot} \cdot ,{d_N})$ is maximized if and only if $|{g({d_1}, \cdot{\cdot} \cdot ,{d_N})} |$ is maximized. Notice that

Apply triangular equation $\cos \alpha - \cos \beta ={-} 2\sin \frac{{\alpha + \beta }}{2}\sin \frac{{\alpha - \beta }}{2}$, then

The equal sign of equation (13) holds if and only if the poling direction of the number m crystal domain is ${d_m} = sign\{ - \sin [\Delta {k_1}\Delta x(m - \frac{1}{2})]\sin (\Delta {k_1}\Delta x/2)\} $, moreover, ${d_m} = sign\{ - \sin [\Delta {k_1}\Delta x(m - \frac{1}{2})]\} $ since $\sin (\Delta {k_1}\Delta x/2) > 0$. Hence the APPLN is degenerated to PPLN with a poling period of $\frac{{2\pi }}{{\Delta {k_1}}}$. It can be understood that since Δk₁=Δk₂ at the degenerate point, a PPLN with identical poling periods of $\frac{{2\pi }}{{\Delta {k_1}}}$ satisfies phase matching of two processes simultaneously.

For understanding the mechanism of CTHG d_eff enhancement at λ₁ = 1420 nm, where Δk₂/Δk₁ = 3, we visualize the local gains of λ₂ and λ₃ 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 ${\Lambda _1}$ is continuously increasing, meanwhile that of SFG also have some increasement after one poling period ${\Lambda _1}$ because of the inversion of nonlinearity after $\frac{{{\Lambda _1}}}{2}$(which is $\frac{{3{\Lambda _2}}}{2}$).

 

Fig. 4. At the fundamental wavelength of 1420 nm, Δ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.

Download Full Size | PDF

The optimal crystal domain distributions, searched by SA algorithm, for Δk₂/Δk₁ = 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 d_eff corresponding to the optimal APPLN crystal near the fundamental wavelength of 1420 nm appears a peak in Fig. 2.

 

Fig. 5. Optimal CTHG crystal with fundamental wavelength of 1420 nm(a), 2080nm(b), and 3550 nm(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 Supplement 1.

Download Full Size | PDF

4. 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, Δk₁ and Δk₂. 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

Fig. 6 shows the comparison of d_eff between the optimal APPLN calculated by Zhang's formula and those by the SA algorithm when the length of the crystal is L = 0.5 mm and 1 mm. 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 λ₁= 1420 nm and λ₁= 3550 nm. In these areas the APPLN superlattice arrangements searched by SA have higher d_eff than those by Zhang's formula. We will discuss these through analysis and simulation in the following.

 

Fig. 6. Comparison of the d_eff 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).

Download Full Size | PDF

The basic idea of our calculation is to extract the sum of the terms related to the crystal domain d_m in the objective function (7), denoted as f_m, and by calculation we show that f_m 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, f_m > 0 is needed to obtain the poling direction of d_m. From the expression of d_m 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

Neglecting the coefficient $\frac{2}{{{L^2}\Delta {k_1}\Delta {k_2}}}$, sum of the terms in (15) associated with d_m is

Since $\cos \alpha - \cos \beta ={-} 2\sin \frac{{\alpha + \beta }}{2}\sin \frac{{\alpha - \beta }}{2}$,

Then (16) can be written as

If the sign of the m^th domain

Comparison of (14) and (20) shows that weight factors of the poling period for SHG and SFG are

In section 3, the optimal APPLN for the fundamental wavelength from 1000 nm to 6000 nm 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 λ₁ = 1000 nm, 3300 nm, 3550 nm, 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.

 

Fig. 7. (a)Weight curves derived reversely from the optimal APPLNs searched by SA at λ1 = 1000 nm, 3300 nm, 3550 nm, respectively.(b)The APPLN domain distributions corresponding to Zhang's analytical formula with L = 0.5 mm. (c) The APPLN domain distribution from the SA optimization.

Download Full Size | PDF

As shown in Fig. 7, the weight ratio curve c(m) at λ₁ = 1000 nm 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 d_eff obtained by (7) is higher than that from Zhang’s formula as shown in Fig. 6(a) at 3300 nm. We think the APPLN designed by us is better than Zhang’s formula at this fundamental wavelength. At degenerate point of 3550 nm, 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, Δk₁=Δk₂, any weight ratio curve gives the result, a PPLN.

It is easy to expect different d_eff 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 3370 nm.

 

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.

Download Full Size | PDF

Actually, the method we proposed above is not only suitable for CTHG (SHG + SFG), but also any coupled two ${\chi ^{(2)}}$ processes (with corresponding Δk₁ and Δk₂) in a single QPM crystal, such as SHG + DFG, DFG + SFG, etc.

5. Upper bond of d_eff

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

Maximizing objective function can be converted to a binary quadratic programming (BQP):

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 d_m in equation (7) only to be +1 and -1, and “relax” d_m into an N-dimensional vector that satisfies |x_m|=1. By applying the semidefinite programming, the original problem can be transformed into

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 ${x_m} = {({q_m},0,0, \cdot{\cdot} \cdot ,0)^T}$, ${q_m} \in \{ + 1, - 1\} $, 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 d_eff comparing to Δx = 0.1µm used in Fig. 2, but benefit to reduce computing time.

 

Fig. 9. (a) Maximum d_eff 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 λ₁ = 1200 nm (b) and λ₁ = 3400 nm (c).

Download Full Size | PDF

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 (λ₁= 3550), i.e. numerical result shows that each x_m corresponding to the solution of problem (23) has only 2 non-zero components, which can be expressed as ${x_m} = {({p_m},{q_m},0, \cdot{\cdot} \cdot ,0)^T}$, where $p_m^2 + q_m^2 = 1$. Fig. 9 (b) illustrates one of the cases, when the fundamental wavelength λ₁ = 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.

6. 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 1000 nm to 6000 nm. The simulation results show that the d_eff curve has peaks at Δk₂/Δk₁= 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 (Δk₂=Δk₁), where the d_eff 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 d_eff 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.